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1.0  SUMMARY 


Relative  orbital  elements  provide  a  geometric  description  of  satellite  motion  referenced  to 
another  satellite  in  close  formation  flight  under  two-body  circular  reference  conditions.  When  the 
relative  motion  is  adequately  characterized  by  first  order  theory  such  as  the  linear  Clohessy- 
Wiltshire  model  and  analytic  solution,  the  six  relative  orbital  elements  describe  the  location  and 
size  of  the  instantaneous  in-the-plane  ellipse,  the  angular  position  around  the  ellipse,  and  the  out- 
of-plane  amplitude  and  angular  position.  These  elements  are  explicitly  relatable  to  the  six 
rectangular  position  and  velocity  coordinates.  Second  order  theories,  such  as  the  Quadratic 
Volterra  model  and  analytic  solution,  become  necessary  when  the  close-proximity  assumption  is 
violated.  Two  distinct  theories  for  (quasi)  second  order  relative  orbital  element  sets  are  explored. 
One  theory  uses  the  expanded  solution  form  and  introduces  several  instantaneous  ellipses  with 
corresponding  element  sets  totaling  twenty-one  individual  elements.  Overall  motion  is  described 
by  a  linear  combination  of  the  ellipses.  Another  theory  uses  the  compacted  solution  form  and 
introduces  a  single  instantaneous  ellipse  with  corresponding  element  set  totaling  nine  individual 
elements.  In  each  case,  the  theory  quantifies  distortion  of  the  first  order  relative  orbital  elements 
when  including  second  order  effects.  The  new  variables  are  described  as  quasi  elements  because 
they  form  a  non-minimal  set  of  coordinates,  and  hence  cannot  be  constructed  as  output  from  a 
general  transformation  having  a  six-dimensional  state  epoch  input,  without  knowledge  of  the 
approximate  analytic  solution.  Several  examples  are  presented,  both  analytically  and 
numerically.  New  periodic  conditions,  geometric  interpretations,  and  relations  to  fixed  and 
variable  Lissajous  curves  and  the  double-pendulum  harmonograph  are  also  explored. 
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2.0  INTRODUCTION 


Although  a  mathematical  representation  of  a  physical  system  can  exhibit  certain  invariant 
behaviors,  many  commonly  encountered  properties  are  dependent  on  the  coordinates  used  to 
represent  the  system.  Such  properties  include  stability,  nonlinearity,  symmetry,  diagonality, 
singularity,  and  others.  Typical  engineering  analysis  efforts  involve  establishing  such  properties 
using  a  given  set  of  coordinates,  followed  almost  involuntarily  (and  incorrectly)  by  replacing  the 
notion  of  coordinate-dependent  properties  with  a  notion  of  system-dependent  properties.  Refs. 
[1,2]  give  a  deep,  enlightening  discussion  of  this  topic,  particularly  for  orbital  mechanics  and 
space  attitude  dynamics.  To  further  complicate  matters,  an  infinite  number  of  coordinate  sets  for 
the  engineer  to  choose  from  exist.  Ref  [2]  develops  a  quantifiable  measure  to  identify  optimum 
coordinate  selection  regarding  the  specific  property  of  nonlinearity,  or  the  minimization  thereof. 
An  astrodynamics  example  investigated  a  coordinate  comparison  of  Classical  Cartesian,  Burdet 
Regularization,  Kustaanheimo-Stiefel  Regularization,  and  Quaternion  Regularization  sets;  while 
a  rotational  motion  example  compared  various  attitude  descriptions  using  Euler  Angle,  Classical 
Rodrigues,  Modified  Rodrigues,  and  Euler  Parameter  sets.  For  the  specific  propagated  cases, 
Kustaanheimo-Stiefe  Regularization  (without  the  temporal  equation)  vs.  Classical  Cartesian 
(with  the  temporal  equation)  and  Euler  Parameter  coordinate  sets  exhibited  the  most  linear 
behavior.  When  performing  coordinate  selections,  further  difficulties  arise  when  emphasis  is 
given  to  the  important  property  of  geometry,  or  the  facilitation  of  insight  thereof  An  analyst  will 
intuitively  recognize  when  one  set  of  coordinates,  through  some  visualization  or  mechanism, 
offers  qualitative  advantages  over  other  coordinate  sets,  yet  a  suitable  procedure  for  making 
choices  based  on  this  property  does  not  appear  to  exist.  This  work  focuses  on  retaining  this  type 
of  geometric  property  for  a  specific  set  of  coordinates  when  increasing  analytic  solution  fidelity 
for  relative  orbital  motion  between  two  satellites.  Throughout  this  discussion,  the  reference 
satellite,  or  reference  point  in  space-time,  is  denoted  "chief,  while  the  perturbed  satellite  is 
called  "deputy". 

For  absolute  motion,  a  large  number  of  coordinate  sets  have  been  developed  with  varying 
attributes,  mostly  defined  with  respect  to  an  inertial  reference  frame  with  origin  coinciding  with 
the  primary  attracting  body.  Ref  [3]  provides  a  thorough  survey  of  at  least  twenty-two  of  these 
sets.  Some  of  these  coordinates  include  classical,  rectangular  (Cartesian),  geographic,  spherical, 
Delaunay,  Poincare,  Kustaanheimo-Stiefel,  equinoctial,  quaternion,  and  many  others  arising 
from  modifications-combinations  of  Euler  angles,  Euler  parameters,  and  functions  of  classical 
elements.  Refs.  [4-6]  provide  other  sets  when  concepts  like  universal  coordinates  and 
fundamental  invariants  are  considered.  Driving  motivations  for  the  creation  of  this  plethora  of 
coordinate  sets  include  singularity  elimination,  dynamical  theory  heritage,  conic  curve 
unification,  and  computational  performance.  Adoption  of  such  coordinates  often  comes  at  the 
expense  of  losing  transparent  geometry,  insight,  and/or  simplicity.  For  relative  motion,  a  smaller 
but  still  large  library  of  coordinate  sets  exists,  which  are  typically  defined  with  respect  to  the 
local-vertical  local-horizontal  reference  frame  with  origin  coinciding  with  the  chief  satellite.  A 
thorough  survey  of  coordinates,  models,  and  solutions  for  relative  motion  is  given  in  Ref  [7]. 
Rectangular,  cylindrical,  and  spherical  coordinate  sets  have  been  considered  at  least  as  far  back 
as  1959-1960,  1961,  1963,  [8-11]  with  lunar-geo  relative  motion  in  a  Cartesian  description  going 
back  to  Hill  [12].  Non-standard  cylindrical  and  spherical  coordinate  definitions  are  commonly 
used  and  based  on  differences  of  absolute  deputy  and  chief  variables  [13].  These  three  sets 
facilitate  strong  connection  to  the  underlying  geometry.  Cylindrical  and  spherical  coordinates 
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possess  the  extra  advantage  of  not  requiring  small  azimuth  angle  assumptions  when  deriving 
simplified  governing  differential  equations  [14,15].  Orbital  element  differenees  (OED)  are  a 
popular  coordinate  choice  and  were  originally  suggested  by  Hill  [12],  These  coordinates  are 
based  on  differencing  absolute  classical  elements  of  the  deputy  and  chief,  but  any  set  of  elements 
can  also  be  employed.  Refs.  [16,17]  first  seriously  considered  this  approach  using  Lagrange 
coefficients  or  F  and  G  functions.  Ref  [18]  used  this  approach  to  analyze  relative  motion  along 
highly  elliptic  orbits.  Refs.  [19,20]  use  OED  based  on  classical  elements  to  generate  higher  order 
Cartesian  solutions  for  elliptic  orbits.  Ref  [21]  formulates  relative  orbit  geometry  and 
relationships  between  Cartesian  coordinates  and  classic  OED  for  first  order  dynamics.  Refs. 
[22,23]  use  OED  based  on  Deprit-Rom  like  elements  [24]  to  remove  singularities  associated  with 
circular  orbits,  and  on  equinoctial  elements  to  remove  singularities  associated  with  both  circular 
orbits  and  equatorial  orbits.  This  formulation  is  known  as  the  geometric  method.  Ref  [25]  uses 
OED  based  on  generalized  relative  eccentricity  and  inclination  vector  separation  coordinates  to 
quantify  deputy  relative  motion,  which  was  used  heavily  in  the  PRISMA  mission.  Ref  [26]  uses 
quaternion-based  elements  to  describe  deputy  motion  that  is  singularity  free,  at  the  expense  of 
increased  complexity  and  indirect  ties  to  the  geometry.  One  of  the  most  intuitive  coordinate  sets 
for  relative  motion  is  relative  orbital  elements  (ROE),  defined  in  Refs.  [27-30].  These 
coordinates  were  conceived  from  the  first  order  circular  chief  analytic  solution  (Clohessy- 
Wiltshire  or  CW  solution)  expressions  for  relative  motion,  and  are  closely  aligned  with  classical 
orbital  elements  in  absolute  motion,  in  that  the  elements  are  explicitly  tied  to  the  solution  path 
geometry.  The  six  relative  orbital  elements  describe  the  location  and  size  of  the  instantaneous  in- 
the-plane  ellipse,  the  angular  position  around  the  ellipse,  and  the  out-of-plane  harmonic 
amplitude  and  angular  position,  and  are  explicitly  relatable  to  the  six  rectangular  position  and 
velocity  coordinates.  Finally,  Ref  [31]  investigated  how  these  relative  orbital  elements  can  be 
extended  under  the  second  order  circular  chief  analytic  solution  (Quadratic  Volterra  or  QV 
solution)  expressions  from  Ref.  [32]. 

The  objective  of  this  work  is  to  further  investigate  extensions  of  relative  orbital  elements 
under  the  QV  solution  framework.  Two  distinct  theories  for  second  order  relative  orbital  element 
sets  are  explored.  One  theory  uses  the  expanded  solution  form  and  introduces  several 
instantaneous  ellipses  with  corresponding  element  sets  totaling  twenty-one  individual  elements. 
Overall  motion  is  described  by  a  linear  combination  of  the  ellipses.  This  theory  is  similar  to  Ref. 
[31]  with  further  development.  Another  theory  uses  the  compacted  solution  form  and  introduces 
a  single  instantaneous  ellipse  with  corresponding  element  set  totaling  nine  individual  elements. 
In  each  case,  the  theory  quantifies  distortion  of  the  first  order  relative  orbital  elements  when 
including  second  order  effects.  Several  examples  are  presented,  both  analytically  and 
numerically.  New  periodic  conditions,  geometric  interpretations,  and  relations  to  fixed  and 
variable  Lissajous  curves  and  the  double-pendulum  harmonograph  are  also  explored. 
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3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 
3.1  FIRST  ORDER  RELATIVE  ORBITAL  ELEMENTS 

Fig.  1  describes  the  relative  motion  of  a  deputy  satellite  with  respect  to  a  chief  satellite.  Two 
reference  frames  are  shown.  The  first  frame  is  the  standard  XYZ  Earth-centered  inertial  (ECI) 
frame.  The  second  frame  is  the  xyz  local-vertical  local-horizontal  (LVLH)  frame  with  origin 
located  on  the  chief,  x  axis  in  the  radial  direction,  y  axis  in  the  transverse  (along-track)  direction, 
and  z  axis  in  the  normal  (cross-track)  direction.  Chief  and  deputy  absolute  position  vectors  in 
algebraic  form  are  denoted  by  Rc  and  Rd  while  the  relative  position  vector  of  deputy  with  respect 
to  chief  is  r.  Position  vector  components  expressed  in  the  LVLH  frame  are 

R,=  [R„OoF 

R  j  =  R,  +  r  =  [(Rc+x)  y  zf 


The  linear  closed-form  CW  relative  trajectory  expressions  for  positions  x(t),  y(t),  and  z(t), 
which  are  accurate  to  first  order,  are  given  in  Eq.  (2)  [9].  In  these  equations,  subscript  "0" 
denotes  the  value  of  a  quantity  at  the  initial  time  to,  and  elapsed  time  from  the  reference  is  x  =  t  - 
to.  Additional  variables  in  Eq.  (2)  include  chief  mean  motion  no,  and  S"  =  sin("),  C"  =  cos("). 

A(t)  ={«Co^,}xo  +  {^o^,}io+{^a^n„t)}yo 

y(t)  =  {6(S  n^T.-n  oT)}xo  +  {l}yo  +  {^-l+CnQ<t)}io  +  DoT-^OoT)}?  o  (2) 

z(0  ={CnoT}zo  +  {^iiox}zo 

Using  trigonometric  identities,  expressions  in  Eq.  (2)  can  be  converted  to  the  form  in  Eq.  (3) 
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(3) 


x(t)  —  Xg(T)  +  — Xg  ^e^a(T) 

y(t)  =  y  c('e)  +  ayS  9y(,)  =  ye(T)  +  2aeS  „(,) 

z(t)  =Zc(-c)  +  agS9^(T)  =beSp(T) _ 

x(t)  =Xc(T)-aje3,(T)S9^(^) 

Kt)  =  +  aydy(T)C0y(,)  =^e('')  +  2aenoCo(,) 

^0  =  Zg(T)  +  a20z(T)C0^T;)  =^e“o^P(x) 

where 

*e  =  4xo  +  ^yo  ,  ye(T)  =  (yo  -  ^xq)  -  (6x9  +  ^yo)noT 

Dq  Hq  “0 

»c  =(3»o+f-yo)^+{i-*o)^ .  “(y)  =noy+‘“'‘«;r*0>'P*0+^<i» 

be  =  (zo)^  +  (^io)^  ,  P(T)  =noT+tan-^{(zo)/(izo)} 

“o  “0 

where  it  becomes  transparent  that  deputy  motion  is  described  by  the  sum  of  an  in-the-plane  yxx 
=  2x1  ellipse  that  drifts  along  the  y  axis  and  a  harmonic  out-of-plane  motion.  Following  Ref. 
[29],  six  ROE  are  defined  directly  from  the  solution  form  in  Eq.  (3)  and  include  constant 
variables  Xe,  ae,  be  and  linear  increasing  variables  ye,  a,  p.  Variables  Xe,  ye  represent  the  in-plane 
ellipse  center  point,  ae  represents  the  in-plane  ellipse  semi-minor  axis,  a  represents  angular 
position  around  the  in-plane  ellipse,  be  represents  the  out-plane  harmonic  amplitude,  and  P 
represents  angular  position  along  the  out-plane  harmonic  path.  These  six  elements  form  a 
minimal  coordinate  realization  of  the  relative  motion.  Eq.  (3)  is  the  transformation  from  ROE  to 
rectangular  states,  while  Eq.  (4)  is  the  reverse  transformation  [29]. 

Xe  =  4x(t)  +  ^y(t)  .  ye(T)  =y(t)  -  :^x(t) 

a?  =  [3x(t)  +  ^y(t)]2  +  [^i(t)]2  ,  a(T)  =tan-l{[^x(t)]  /  [3x(t)  +  ^y(t)]}  (4) 

Dq  Hq  IIq  Qq 

be  =  [z(t)]^  +  [^2(t)]^  ,  P(T)  =tan-l{[z(t)] /  [^i(t)I} 

“o  “o 
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3.2  FIXED  AND  VARIABLE  LISSAJOUS  CURVES 


Before  introducing  second  order  ROE,  a  discussion  of  Lissajous  curve  theory  and  the  double¬ 
pendulum  harmonograph  is  presented.  These  topics  are  closely  related  to  the  concept  of  second 
order  ROE.  A  two-dimensional  Lissajous  curve  is  generated  by  the  parametric  equations  below 
[33-35]. 


x(t)  =aj^cos(tOx'c-(|)x) 
y(t)  =aySin(a)yX-wl)y) 

These  curves  were  first  studied  by  Bowditch  in  1815  and  later  by  Lissajous  in  1857.  Eq.  (5)  can 
represent  a  wide  variety  of  coupled  harmonic  responses  to  initial  excitation  or  harmonic  forcing 
in  physical  systems,  and  the  x(t)  vs.  y(t)  Lissajous  curves  they  produce  are  commonly 
encountered  during  experimental  signal  analysis  with  modem  instmmentation.  Elementary 
curves  can  be  generated  by  certain  combinations  of  parameters:  cox  =  coy  and  (|)x  +  (|)y  =  ±7r/2 
produces  a  line;  cox  =  2coy  and  (jix  +  2^y  =  0,  ±71,  or  2cox  =  coy  and  2(|)x  +  (|)y  =  ±7i/2,  produces  a 
parabola;  cox  =  coy,  ([ix  +  ([ly  =  0,  ±71  and  ax  =  ±ay  produces  a  circle;  cox  =  coy,  (|)x  +  (|)y  =  0,  ±71  and  ax 
±ay  produces  an  unrotated  ellipse;  and  cox  =  tOy,  (|)x  +  (|)y  0,  ±7i:/2,  ±71  and  ax  7^  ±ay  or  ax  =  ±ay 

produces  a  rotated  ellipse.  Although  these  simple  cases  are  examples  of  Lissajous  curves,  the 
implied  meaning  of  Lissajous  curves  is  often  reserved  for  more  complex  shapes.  Multi-lobe 
curves,  or  equivalently  multiple  extrema  points,  are  generated  when  frequency  ratio  cox/coy 
consists  of  higher  integer  values,  such  as  3/1  or  3/2.  Finally,  closed  Lissajous  curves  occur  when 
(Ox/coy  is  rational  (commensurate),  and  when  irrational  (non-commensurate),  the  curves  are  open. 

Suppose  a  particle  is  "orbiting"  the  center  point  of  a  rotated  elliptic  Lissajous  path  as 
indicated  in  Fig.  2  for  the  specific  set  of  parameter  values.  In  general,  four  state  variables  are 
needed  to  define  the  planar  particle  motion.  Along  this  programmed  path,  four  suitable  states  or 
"orbital"  elements  are  ax,  ay,  0x(t)  =  coxX  -  (|)x,  0y(T)  =  coyX  +  (jiy;  the  first  two  elements  describe  the 
path  size-shape,  the  last  two  elements  describe  the  particle  position.  Note  this  situation  is  more 
involved  than  the  first  order  CW  curvilinear  motion  component  around  the  instantaneous  ellipse 
since  no  single  underlying  parameter  exists  to  determine  the  motion  amplitudes  and  also  there  is 
no  common  angle  for  the  harmonic  function  arguments.  To  further  underscore  how  ax,  ay,  0x(t), 
0y(T)  can  be  used  for  motion  interpretation,  these  elements  are  shown  in  Fig.  2  along  with  two 
auxiliary  circles  of  radius  ax  and  ay.  Fixed  elements  ax  and  ay  describe  dimensions  of  a 
rectangular  area  that  encloses  the  orbit  path.  Note  the  path  is  tangent  to  the  boundary  of  the 
rectangular  area  at  four  points,  not  generally  coincident  with  the  principal  axis  points  of  the 
rotated  ellipse.  The  particle  initial  position  at  time  t  =  to  or  x  =  0,  and  a  general  position  at  time  t 
or  X,  are  highlighted  in  Fig.  2.  The  actual  particle  on  the  Lissajous  path  projects  to  shadow 
particles  moving  along  the  auxiliary  circles,  also  indicated  in  Fig.  2.  Linearly  increasing  elements 
0x(x)  and  0y(x)  describe  the  angular  position,  measured  from  the  x  axis,  of  the  shadow  particles 
around  the  auxiliary  circles  of  constant  radii  equal  to  elements  ax  and  ay,  which  in  turn  locates  the 
actual  particle  on  the  Lissajous  path.  Another  but  equivalent  viewpoint  is  where  angles  ooxx  and 
coyX  locate  the  shadow  particles  from  the  reference  marks  determined  by  phase  angles  (t)x,  (t)y, 
which  correspond  to  the  initial  position  at  t  =  to  or  x  =  0.  The  four  elements  and  the  auxiliary 
circles  through  Eq.  (5)  provide  a  clear  depiction  of  the  particle  motion.  With  this  geometry  and 
Eq.  (5),  a  clear  representation  of  the  motion  is  captured.  Given  ax  and  0x(x)  (shadow  particle  on 
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the  X  auxiliary  circle),  the  actual  particle  and  its  x  coordinate  can  always  be  generated  and 
interpreted  by  a  cosine  projection.  Likewise,  given  ay  and  0y(T)  (shadow  particle  on  y  auxiliary 
circle),  the  actual  particle  and  its  y  coordinate  can  always  be  generated  and  interpreted  by  a  sine 
projection. 


Figure  2.  "Orbiting"  Particle  on  Fixed  Lissajous  Curve  (In-Plane  View) 

The  two-dimensional  Lissajous  curve  considered  thus  far  can  be  classified  as  fixed,  meaning 
ax,  ay,  cox,  ooy,  (|)x,  (|)y  arc  constant.  A  variable  Lissajous  curve  is  now  defined  where  these 
parameters  are  time  dependent.  A  two-dimensional  variable  Lissajous  curve  is  thus  generated  by 
the  parametric  equations  below. 


x(t)  =ax(x)cos{ci)x(T)x-<l)x(x)} 

y(t)  =ay(x)sin{(i>y(x)x+(l>y(x)}  ^  ^ 

Variable  Lissajous  curves  provide  a  considerably  more  diverse  framework  that  is  well-suited  for 
characterizing  higher  order  approximations  of  relative  orbital  motion.  The  concept  of  a  variable 
Lissajous  curve  is  not  necessarily  new.  Fig.  3  shows  the  double-pendulum  harmonograph,  a 
mechanical  apparatus  consisting  of  two  separated  pendulums  swinging  about  orthogonal  axes, 
with  an  attached  marker  that  records  on  paper  a  complex  curve  [36].  This  device  could  be 
considered  an  analog  propagator  for  physical  systems  that  exhibit  multi-axis  complex  harmonic 
motion,  including  relative  satellite  motion.  The  concept  can  also  be  mechanized  by  harmonic 
signals  driving  a  gimbal  mounted  laser  to  trace  out  the  curve.  The  harmonograph  is  quite  popular 
in  mathematical  description  of  tones,  music,  and  sound,  along  with  generating  modem  artistic 
drawings  or  visual  entertainment  [37,38].  Under  ideal  assumptions  with  no  pivot  friction  or  air 
resistance,  small-amplitude  pendulum  motion,  and  a  stiff  pivot  support,  the  harmonograph 
generated  figure  is  a  fixed  Lissajous  curve  governed  by  Eq.  (5).  An  example  of  a  variable 
Lissajous  curve  is  generated  when  non-ideal  assumptions  are  considered:  frictional  damping 
leading  to  variable  amplitude,  large-angle  swinging  leading  to  variable  frequency,  and  soft 
supporting  leading  to  variable  phase.  Various  sources  describe  the  case  where  ax(T)  and  ay(T)  are 
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exponential  decaying  functions  reflecting  system  damping.  With  servo  driven  pendulums,  the 
harmonograph  could  be  made  to  generate  tailored  multi-axis  complex  harmonic  motions. 


Figure  3.  Harmonograph  with  Generated  Curve 

(Figure  courtesy  of  Paul  Bourke,  http://paulbourke.net/geometry/harmonograph/) 


Now  suppose  the  particle  is  "orbiting"  the  center  point  of  the  variable  Lissajous  path,  a 
continuously  varying  and  distorting  rotated  planar  ellipse.  In  general,  four  state  variables  are 
again  necessary  and  sufficient  to  define  the  planar  particle  motion.  Four  suitable  "orbital" 
elements  are  ax(T),  ay(T),  0x(t)  =  cox(t)t  -  (|)x(t),  0y(T)  =  coy(T)T  +  (t)y(T);  which  have  the  same 
geometric  meaning  as  before,  except  that  amplitude,  frequency,  and  phase  are  varying  with  time. 
Fig.  4  illustrates  this  geometry.  The  actual  particle  on  the  variable  Lissajous  path  projects  to 
shadow  particles  moving  along  the  variable  radii  instantaneous  auxiliary  circles  (ax(T),ay(T))  with 
variable  rates  (oc)x(T),oc)y(T))  measured  from  variable  reference  marks  ((|)x(T),(|)y(T)).  This  concept  is 
similar  to  the  concept  already  employed  with  first  order  ROE  where  the  relative  orbital  path  is 
described  by  a  continuously  moving  or  instantaneous  ellipse.  Here,  the  conceptual  elements  are 
instantaneous  auxiliary  circles  and  the  instantaneous  "fixed"  Lissajous  curve.  Given  ax(T)  and 
0x(t)  (shadow  particle  on  the  x  auxiliary  "circle"),  the  actual  particle  and  its  x  coordinate  can 
always  be  generated  and  interpreted  by  a  cosine  projection.  Likewise,  given  ay(T)  and  0y(T) 
(shadow  particle  on  y  auxiliary  "circle"),  the  actual  particle  and  its  y  coordinate  can  always  be 
generated  and  interpreted  by  a  sine  projection.  With  the  Fig.  4  geometry  and  Eq.  (6),  a  clear,  but 
more  intellectually  involved,  representation  of  the  motion  is  captured.  All  four  time  varying 
elements  are  physical  variables  that  can  be  indicated  on  a  drawing  of  the  motion,  and  hence  are 
geometrically  meaningful.  This  framework  is  also  very  similar  to  the  idea  of  osculating  orbital 
elements  in  absolute  satellite  motion  used  to  preserve  the  geometric  usefulness  of  the  coordinates 
when  higher  fidelity  modeling  is  incorporated,  i.e.,  perturbations. 
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Figure  4.  "Orbiting"  Particle  on  Variable  Lissajous  Curve 


Several  examples  of  variable  Lissajous  curves  are  studied  next  to  illustrate  some  of  the 
features  and  signatures  these  curves  can  exhibit.  Only  variable  phase  and  variable  amplitude 
cases  are  considered,  since  variable  frequency  relative  orbital  motion  does  not  seem  to  occur,  at 
least  for  the  assumptions  commonly  used  in  deriving  approximate  higher  order  analytic  solutions 
[32],  Using  the  parameter  values  from  Fig.  2  as  a  baseline.  Fig.  5  shows  an  example  variable 
phase  Lissajous  curve  with  a  linearly  time  varying  phase  function  (|)x(t).  Phase  angle  (|)x(t)  starts 
at  +71/8  and  decreases  by  7i/4  for  every  cycle  of  motion.  Four  cases  are  shown  with  final  times  of 
1(27i),  2(271),  4(27i),  8(271).  For  a  final  time  of  l(27i),  the  particle  returns  to  the  initial  position 
forming  an  unsymmetric  closed  curve.  Closed  asymmetry  occurs  when  2cos(-7i/8)  =  x(0)  = 
x(  1(271:))  =  2cos(+7r/8).  For  a  final  time  of  2(27i:),  note  how  the  continuously  varying  ellipse  is 
rotating  clockwise  as  (t)x(T)  decreases,  or  -(|)x(t)  increases.  In  the  context  of  Fig.  4,  the  reference 
line  corresponding  to  -(|)x(t)  is  advancing,  causing  the  x  harmonic  to  peak  at  an  earlier  time  for 
each  cycle,  which  has  the  effect  of  rotating  the  variable  ellipse  in  a  retrograde  sense.  If  an 
auxiliary  circle  was  introduced,  the  shadow  particle  would  be  located  for  this  changing  angle 
through  element  0x(t)  =  coxx  -  (|)x(t),  and  the  x  coordinate  of  the  actual  particle  could  be  easily 
determined.  For  a  final  time  of  8(27i),  the  particle  returns  to  the  initial  position  forming  an 
asymmetric  closed  curve:  2cos(-7i/8)  =  x(0)  =  x(8(27i))  =  2cos(8(27r)+157r/8)  =  2cos(157i/8)  = 
2cos(-7r/8). 

Again  using  the  parameter  values  from  Fig.  2  as  a  baseline.  Fig.  6  shows  an  example  variable 
amplitude  Lissajous  curve  with  a  linearly  time  varying  amplitude  function  ax(x).  Amplitude  ax(T) 
starts  at  2  and  decreases  by  0.2  for  every  cycle  of  motion.  Four  cases  are  shown  with  final  times 
of  1(271:),  2(271),  4(2 7i),  8(27r).  For  a  final  time  of  l(27i:),  the  particle  closely  approaches  but  does 
not  reach  the  initial  position  forming  an  open  curve.  Openness  occurs  because  2cos(-7i/8)  =  x(0) 
x(1(27i))  =  L8cos(-7i/8).  For  a  final  time  of  2(27i),  note  how  the  continuously  varying  ellipse  is 
shrinking  in  width  as  ax(T)  decreases.  Referring  back  to  Fig.  4,  the  auxiliary  circle  radius  related 
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Figure  5.  Variable  Phase  Lissajous  Curve 

to  ax(T)  is  decreasing,  causing  the  x  harmonic  to  diminish  over  each  cycle  that  occurs,  which  has 
the  effect  of  initially  circularizing  the  variable  ellipse  (final  time  of  4(27i))  and  then  leading  to  a 
reversal  of  the  semi-major  and  semi-minor  axes  (final  time  of  8(27i:)).  If  an  auxiliary  circle  was 
introduced,  the  shadow  particle  would  be  located  for  this  changing  radius  through  element  ax(T), 
and  the  x  coordinate  of  the  actual  particle  could  be  easily  determined.  In  this  variable  amplitude 
case,  the  Lissajous  curve  never  becomes  a  closed  curve. 

As  a  final  two-dimensional  example.  Fig.  7  shows  an  example  variable  phase  and  variable 
amplitude  Lissajous  curve  with  a  linearly  time  varying  phase  function  (|)x(t)  and  linearly  time 
varying  amplitude  functions  ax(T)  and  ay(T).  Phase  angle  (|)x(t)  starts  at  and  decreases  by 
71/16  for  every  cycle  of  motion.  Amplitude  ax(t)  starts  at  2  and  increases  by  0.2  for  every  cycle  of 
motion,  while  amplitude  ay(T)  starts  at  1  and  increases  by  0.1  for  every  cycle  of  motion.  Four 
cases  are  shown  with  final  times  of  l(27i),  2(27i:),  4(27i),  8(27i).  For  a  final  time  of  l(27i:),  the 
particle  approaches  but  does  not  reach  the  initial  position  forming  an  opened  curve.  Openness 
occurs  because  2cos(-7r/8)  =  x(0)  x(l(27i))  =  2.2cos(-7i/16)  and  lsin(-i-7r/8)  =  y(0)  y(l(27i))  = 

l.lsin(-i-7r/8).  For  a  final  time  of  2(27i),  note  how  the  continuously  varying  ellipse  is  rotating 
clockwise  as  (|)x(t)  decreases,  or  -(|)x(t)  increases.  In  the  context  of  Fig.  4,  the  reference  line 
corresponding  to  -(t)x(T)  is  again  advancing,  causing  the  x  harmonic  to  peak  at  an  earlier  time  for 
each  cycle,  which  has  the  effect  of  rotating  the  variable  ellipse  in  a  retrograde  sense.  Also  note 
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Figure  6.  Variable  Amplitude  Lissajous  Curve 

the  continuously  varying  ellipse  is  expanding  as  ax(T)  and  ay(T)  increase  in  a  2-to-l  fashion.  In 
the  context  of  Fig.  4,  the  auxiliary  circle  radii  corresponding  to  ax(T)  and  ay(T)  are  increasing, 
causing  the  x  and  y  harmonics  to  swell  over  each  cycle  that  occurs,  which  has  the  effect  of 
preserving  the  2x1  aspect  ratio  as  the  variable  ellipse  rotates  (final  times  of  4(27i)  and  8(27r)).  If 
an  auxiliary  circle  was  introduced,  the  shadow  particles  would  be  located  for  these  changing 
distances  and  angles  through  elements  ax(T),  ay(T),  0x(t)  =  coxX  -  (|)x(t),  0y(T)  =  coyX  +  {|)y,  and  the 
x,y  coordinates  of  the  actual  particle  could  be  easily  determined.  In  this  example,  the  variable 
Lissajous  curve  is  open  and  never  repeats  a  periodic  path. 

A  three-dimensional  fixed  Lissajous  space  curve  is  generated  by  the  parametric  equations 
below  [33,39-41]. 


x(t)  =axCOs((OxT-<|)x) 

y(t)  =aySin(MyT-i-(j)y)  (7) 

z(t)  =a2sin([02T-i-(|)2) 

General  space  curves  can  be  classified  according  to  whether  they  are  closed  vs.  open  (periodic 
vs.  non-periodic)  and  self-intersecting  vs.  non-self-intersecting  (conjunctive  vs.  non¬ 
conjunctive).  Fig.  8  illustrates  the  four  possibilities  and  indicates  required  conditions  for  each 
case. 
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Figure  7.  Variable  Phase- Amplitude  Lissajous  Curve 

From  these  conditions,  closed  Lissajous  space  curves  occur  when  cox/coy,  ooy/ooz,  coz/cox  are  all 
rational,  and  when  at  least  one  of  the  frequency  ratios  is  irrational,  the  curves  are  open.  Further, 
several  statements  can  be  proved  regarding  Eq.  (7)  Lissajous  space  curve  intersections,  all  of 
which  require  cox,  coy,  coz  to  be  pairwise  relative  prime  (greatest  common  denominator  is  1).  First, 

if  (|)y/COy  "t"  (|)x/COx  ^  TT  {(2ny“t“l  j/cOy  -  Ox/OOx},  (|)z/00z  “  (|)y/COy  ^  71  {(2nz"t”l)/C0z  “  (2ny“t“l  j/oOy} ,  -  (|)x/00x  - 

(|)z/coz  7i{nx/cox  -  (2nz+l)/coz}  for  integer  nx,  ny,  nz,  then  the  curve  is  non-self-intersecting. 
Second,  if  no  more  than  two  of  the  previous  phase-frequency  conditions  are  equalities  rather  than 
inequalities,  then  the  curve  is  self-intersecting  with  a  finite  number  of  crossings.  Note  this 
statement  excludes  the  exceptional  case  where  two  of  the  three  frequencies  are  equal  to  1  (which 
still  satisfies  the  relative  prime  condition).  Third,  if  all  three  of  the  previous  phase-frequency 
conditions  are  equalities  rather  than  inequalities,  then  the  curve  is  self-intersecting  with  an 
infinite  number  of  crossings.  Note  this  case  corresponds  to  a  singular  curve  where  the  curve 
backtracks  along  the  same  path  causing  every  point  to  be  a  crossing  point,  although  when  plotted 
the  curve  gives  the  "appearance"  of  a  non-self-intersecting  curve.  A  non-conjunctive  periodic 
Lissajous  curve  is  called  a  Lissajous  knot  in  mathematical  knot  theory  [41].  A  three-dimensional 
periodic  orbit  in  the  restricted  three-body  problem,  although  not  specifically  generated  by  Eq. 
(7),  is  similar  to  a  Lissajous  knot  and  is  commonly  referred  to  as  a  halo  orbit  [42-44].  Further,  a 
three-dimensional  quasi-periodic  (nearly  periodic)  orbit  solving  the  restricted  three-body 
problem  is  referred  to  as  a  Lissajous  orbit  [42-44]. 
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For  specific  t: 

x(0  =  &  y<t)  =  y(t+Al)  &  ^(t)  =  z(t+At) 

for  only  At  =  nT  .  ii  =  ± I  ^2^3  . . . 


For  specific  L,t'  where  t  sa  t': 
xCt)  =  x(t+Al>  &:  y(t)  =  y(t+At)  Sl  z(t)  =  z<t+Al> 
for  only  At  =  nT  ,  n  =  ±  I  ^2^3  , . . 
uud 

x<t')  =  x(t'+At)  &  y(t')  =  y(t'+At3  &  z(t')  =  z(t'+Al) 
for  some  0  <  At  <  T  &  -T  <  At  c  0 


For  speeilie  t, 

x(t>  ^  x<l+At)  &  y(t)  ^  y(l+At)  &  z(t)  ^  z<t+Al> 
for  all  At  >  0  &  At  <  0 


Fiir  specific  tjt'  where  t  sa  t': 
x(  t)  ^  x(t+At)  &  y(t)  y(t+At)  &  z(t>  ztt+At> 
for  all  At  >  0  At  <  0  other  than  At  =  t*  -  t 
and 

x(t')  =  xU’+At)  &  y(t')  =  y(l'+At)  &  z<t')  =  z(t'+At) 
for  some  At  >  0  &  At  <  0  but  not  At  =  nT 


Figure  8.  Space  Curve  Classification 


Figure  9.  "Orbiting"  Particle  on  Fixed  Lissajous  Curve  (Out-Plane  View) 

Suppose  a  particle  is  "orbiting"  the  center  point  of  a  rotated  elliptic  Lissajous  space  path.  The 
projected  in-plane  x,y  motion  appears  as  in  Fig.  2,  and  simultaneously  the  projected  out-plane  x,z 
motion  is  indicated  in  Fig.  9  for  the  specific  set  of  parameter  values.  In  general,  six  state 
variables  are  needed  to  define  the  space  particle  motion.  Along  this  programmed  path,  six 
suitable  states  or  "orbital"  elements  are  ax,  ay,  az,  0x(t)  =  coxX  -  (|)x,  0y(T)  =  coyX  +  {|)y,  0z(t)  =  cozX  + 
{|)z;  the  first  three  elements  describe  the  path  size-shape,  the  last  three  elements  describe  the 
particle  position.  The  x  and  y  related  orbital  elements  can  be  used  for  in-plane  motion  tracking  as 
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previously  discussed,  and  the  x  and  z  related  orbital  elements  can  be  used  in  the  same  way  to 
describe  the  out-plane  motion.  Fig.  9  shows  the  projected  Lissajous  planar  curve  and  auxiliary 
circles  for  the  xz  plane.  Fixed  elements  ax  and  az  describe  dimensions  of  a  rectangular  area  that 
encloses  the  orbit  path.  The  particle  initial  and  general  positions  on  the  Lissajous  path  are 
highlighted  in  Fig.  9,  and  the  corresponding  shadow  particles  appear  on  the  auxiliary  circles. 
Linearly  increasing  elements  0x(t)  and  0z(t)  describe  the  angular  position  of  the  shadow  particles 
around  the  auxiliary  circles,  which  in  turn  locates  via  projections  the  actual  particle  on  the 
Lissajous  path.  The  six  elements  and  the  auxiliary  circles  in  Figs.  2  and  9  along  with  Eq.  (7) 
provide  a  clear  depiction  of  the  particle  motion.  With  this  geometry  and  Eq.  (7),  a  clear 
representation  of  the  motion  is  captured.  Rather  than  interpreting  the  motion  geometry  in  two 
separate  planar  images,  a  three-dimensional  image  can  be  constructed  where  all  three  projections 
are  processed  in  one  view.  However,  for  purposes  of  clarity,  two  separate  three-dimensional 
images,  one  containing  only  the  in-plane  element  geometry  and  the  other  containing  only  the 
out-plane  geometry,  are  shown  in  Fig.  10.  To  further  reduce  clutter  in  Fig.  10,  only  the  particle 
general  position  is  shown  but  the  auxiliary  circles  are  retained.  Computer  generated  three- 
dimensional  color  images  can  facilitate  the  practical  use  of  the  six  "orbital"  elements  and  the 
geometrical  information  they  depict.  In  Fig.  10,  0a  denotes  azimuth  view  angle  measured  from  -y 
axis  with  a  right-hand  rotation  sense  and  0e  denotes  elevation  view  angle  positive  above  xy 
plane. 

The  final  topic  needed  before  returning  to  second  order  ROE  for  deputy-chief  satellite 
motion,  is  variable  Lissajous  space  curves.  A  three-dimensional  variable  Lissajous  curve  with 
time  dependent  amplitude,  frequency,  and  phase  parameters  is  generated  by  the  equation  set 
below. 


x(t)  =ax(x)cos{tOx(T)T-(|)x(x)} 

y(t)  =ay(T)sin{(Oy(T)T+(t)y(T)}  (8) 

z(t)  =a2(x)sin{a32(x)x+(J)2(x)} 

Suppose  the  particle  is  "orbiting"  the  center  point  of  a  variable  Lissajous  space  path,  a 
continuously  varying  and  distorting  rotated  space  ellipse.  In  general,  six  state  variables  are  again 
necessary  and  sufficient  to  define  the  space  particle  motion.  Six  suitable  "orbital"  elements  are 

ax(T),  ay(T),  az(T),  0x(t)  =  cox(t)t  -  (|)x(t),  0y(T)  =  coy(T)T  +  (|)y(T),  0z(t)  =  coz(t)t  +  (|)z(t);  which  have 

the  same  geometric  meaning  as  before,  except  that  amplitude,  frequency,  and  phase  are  varying 
with  time.  Fig.  1 1  illustrates  this  geometry.  The  actual  particle  on  the  variable  Lissajous  path 
projects  to  shadow  particles  moving  along  the  variable  radius  auxiliary  circles  (ax(T),ay(T)  and 
ax(T),az(T))  with  variable  rates  (cox(T),(j0y(T)  and  cox(T),(jOz(T))  measured  from  variable  reference 
marks  ((|)x(T),(|)y(T)  and  (t)x(T),(t)z(T))  in  the  xy  and  xz  planes.  Given  ax(T)  and  0x(t)  (shadow  particle 
on  the  instantaneous  x  auxiliary  circle),  the  actual  particle  and  its  x  coordinate  can  always  be 
generated  and  interpreted  by  a  cosine  projection.  Likewise,  given  ay(T)  and  0y(T)  (shadow 
particle  on  instantaneous  y  auxiliary  circle),  and  az(T)  and  0z(t)  (shadow  particle  on  the 
instantaneous  z  auxiliary  circle),  the  actual  particle  and  its  y  and  z  coordinates  can  always  be 
generated  and  interpreted  by  sine  projections.  In  Fig.  II,  the  instantaneous  "fixed"  Lissajous 
curve  is  a  planar  curve  while  the  varying  Lissajous  curve  is  a  non-planar  twisted  curve  [45]. 
With  the  Fig.  1 1  geometry  and  Eq.  (8),  a  systematic  and  geometric  characterization  of  the  motion 
is  captured. 
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Figure  10.  "Orbiting"  Particle  on  Fixed  Lissajous  Curve  (Perspective  View) 


Figure  11.  "Orbiting"  Particle  on  Variable  Lissajous  Space  Curve 
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Only  one  example  of  a  variable  phase  and  variable  amplitude  Lissajous  space  curve  is 
considered  here.  The  planar  example  in  Fig.  7  serves  as  the  baseline,  and  is  extended  with  the  z 
dimension  motion.  Fig.  12  shows  the  example  curve  with  a  linearly  time  varying  phase  function 
(|)x(t)  and  linearly  time  varying  amplitude  functions  ax(T),  ay(T),  and  az(T).  As  in  the  previous 
example,  phase  angle  (|)x(t)  starts  at  +71/8  and  decreases  by  71/ 16  for  every  cycle  of  motion,  and 
amplitude  ax(T)  starts  at  2  and  increases  by  0.2  for  every  cycle  of  motion,  while  amplitude  ay(T) 
starts  at  1  and  increases  by  0.1  for  every  cycle  of  motion.  New  amplitude  az(T)  starts  at  1  and 
increases  by  0.1  for  every  cycle  of  motion,  which  is  identical  to  ay(T).  Also  the  new  z  axis 
frequency  and  phase  are  constant  and  relate  to  corresponding  y  axis  values  as  ooz  =  coy  and  (|)z  =  - 
(t)y.  Only  one  curve  with  a  final  time  of  8(27i)  is  examined,  but  four  different  perspective  views 
are  shown  in  Fig.  12.  First,  if  the  amplitudes  ax(T),  ay(T),  az(T)  and  phase  (|)x(t)  were  maintained  at 
their  initial  values,  the  Lissajous  curve  would  be  a  closed  curve  (an  ellipse)  lying  in  a  fixed  plane 
skewed  from  the  horizontal,  similar  to  the  curve  shown  in  Fig.  10.  Second,  if  the  amplitudes 
ax(T),  ay(T),  az(T)  varied  as  functionally  indicated  but  phase  (|)x(t)  was  still  maintained  at  the  initial 
value,  the  Lissajous  curve  would  be  an  open  curve  but  still  lying  in  the  fixed  skewed  plane.  The 
curve  would  appear  as  an  expanding  spiral  preserving  the  initial  xxyxz  =  2x1x1  elliptic  shape. 
Finally,  when  the  amplitudes  ax(T),  ay(T),  az(T)  and  phase  (|)x(t)  vary  as  functionally  indicated,  the 
Lissajous  curve  becomes  an  open,  non-self-intersecting,  space  curve,  as  shown  in  Fig.  12.  In  this 
example,  the  twisted  path  characteristic  is  due  to  the  decreasing  phase  (t)x(T).  Expanding  and 
twisting  signatures  are  observed  in  Fig.  12.  In  the  0a  =  30  deg  view,  note  the  growing  positive  z 
and  negative  x  motions  in  the  upper  left  and  the  growing  negative  z  and  positive  x  motions  in  the 
lower  right.  With  this  azimuth  angle,  growth  in  the  y  motion  is  obscured.  Simultaneously,  note 
how  the  peak  x  motions  occur  earlier  for  each  cycle  of  motion  as  the  curve  is  twisted  in  a 
retrograde  sense.  Although  each  image  is  of  the  same  curve,  individual  azimuth  perspectives 
give  different  characteristic  "shapes":  0a  =  30  deg  produces  a  conchiglie  pasta  or  seashell  shape, 
0a  =  +120  deg  produces  a  deformable  membrane  shape,  0a  =  -60  deg  produces  a  cornucopia 
shape,  and  0a  =  -100  deg  produces  a  multi-coil  spring  shape.  The  shape  of  the  Lissajous  curve  is 
really  none  of  these,  and  perhaps  oversimplifying,  the  Lissajous  curve  is  simply  a  complex  three- 
dimensional  curve.  However,  the  concept  of  a  continuously  varying  and  distorting  elliptic 
Lissajous  curve  is  naturally  "visualized"  when  examining  Fig.  12.  Further,  the  geometry  of  this 
osculating  construct  is  precisely  described  by  the  methods  previously  discussed. 
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+30  deg 
+20  deg 


Figure  12.  Variable  Phase-Amplitude  Lissajous  Space  Curve 
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4.0  RESULTS  AND  DISCUSSION 

4.1  SECOND  ORDER  RELATIVE  ORBITAL  ELEMENTS  -  EXPANDED  FORM 

Return  to  the  topic  of  second  order  ROE  for  deputy-chief  motion.  The  nonlinear  closed-form 
QV  relative  trajectory  expressions  for  x(t),  y(t),  and  z(t),  which  are  accurate  to  second  order,  are 
given  in  Eq.  (9)  [32],  In  these  equations,  subscript  "0"  again  denotes  the  value  of  a  quantity  at  the 
initial  time  to,  and  elapsed  time  from  the  reference  is  x  =  t  -  to.  Additional  variables  in  Eq.  (9) 
include  chief  mean  motion  no,  chief  mean  radius  l^o  =  Rc,  and  S"  =  sin("),  C"  =  cos(").  By 
grouping  terms  according  to  individual  expansion  functions:  1,  noX,  (nox)^,  cos(nox),  sin(nox), 
noxcos(nox),  noxsin(nox),  cos(2nox),  sin(2nox),  and  using  trigonometric  identities,  expressions  in 
Eq.  (9)  can  be  converted  to  the  form  in  Eq.  (10),  where  Eq.  (11)  provides  lower  level  coefficient 
definitions.  Eq.  (10)  indicates  trajectory  expressions  contain  three  sets  of  three-dimensional 
parametric  Lissajous  functions:  one  for  harmonic  noX  indicated  by  subscript  "hi",  one  for  the 
product  of  secular  Uox  with  harmonic  noX  indicated  by  subscript  "shl",  and  one  for  harmonic  2noX 
indicated  by  subscript  "h2",  in  addition  to  a  set  of  polynomial  functions  of  noX  indicated  by 
subscript  "c".  The  Lissajous  functions  produce  three  elliptical  space  curves  and  correspondingly 
three  sets  of  xy  and  xz  planar  ellipses,  and  the  polynomial  functions  produce  a  secular  quadratic 
space  curve.  Following  Sections  3-4  and  Ref  [31],  a  set  of  twenty-one  ROE  are  defined  directly 
from  the  expanded  solution  form  in  Eq.  (10)  and  include:  {xc,  yc,  Zc},  {ax, hi,  ay, hi,  az,hi,  0x,hi,  0y,hi, 
0z,hi},  {ax,shi,  ay,shi,  az,shi,  0x,shi,  0y,shi,  0z,shi},  {ax,h2,  ay,h2,  az,h2,  0x,h2,  0y,h2,  0z,h2}.  Elements  {Xc,  yc, 
Zc}  represent  the  in-plane  and  out-plane  center  point  coordinates  for  the  two-dimensional 
ellipses.  Elements  {ax,hi,  ay,hi,  az,hi,  0x,hi,  0y,hi,  0z,hi}  characterize  the  size-shape  of,  and  angular 
position  along,  the  two-dimensional  ellipses  corresponding  to  harmonic  noX.  Recall  each  ellipse 
will  have  two  associated  auxiliary  circles  that  are  explicitly  described  by  the  elements,  which  in 
turn  implicitly  describe  the  actual  geometry  of  the  motion  along  the  ellipse  in  a  meaningful  way. 
Elements  {ax,shi,  ay,shi,  az,shi,  0x,shi,  0y,shi,  0z,shi}  and  {ax,h2,  ay,h2,  az,h2,  0x,h2,  0y,h2,  0z,h2}  play 
similar  roles  for  the  (secular)x(harmonic)  noX  and  harmonic  2noX  ellipses.  Note  the  pure 
harmonic  ellipses  correspond  to  fixed  Lissajous  curves,  but  the  mixed  secular-harmonic  ellipses 
are  linearly  expanding  variable  Lissajous  curves. 

A  comparison  of  second  order  and  first  order  ROE  expressions  is  considered  next.  The  center 
point  coordinate  Xe  is  constant  in  the  first  order  theory  (Eq.  (3))  but  Xc  is  time  dependent  in  the 
second  order  theory  (Eq.  (10)).  This  time  dependency  involves  secular  terms  of  powers  0,1,2 
indicated  by  subscripts  "s0","sl","s2"  in  Eq.  (10).  Observe  the  first  two  terms  in  coefficient  bx,so, 
which  are  linear  functions  of  initial  conditions  xo  and  Yo  (Eq.  (11)),  are  precisely  the  first  order 
Xe  expression  terms.  Other  second  order  initial  condition  terms  in  bx,s0  represent  a  constant 
distortion  to  the  first  order  ROE  Xe.  Other  distortions  that  vary  with  time  are  represented  by  the 

bx, si  and  bx,s2  coefficients  in  Eq.  (10).  In  a  similar  fashion,  the  center  point  coordinate  ye  involves 
static  and  linearly  increasing  terms  in  the  first  order  theory  (Eq.  (3)),  and  this  behavior  for  yc  is 
preserved  in  the  second  order  theory  but  with  distortion  of  the  relevant  parameters.  The  first  two 
terms  in  coefficient  by,s0  involve  initial  conditions  yo  and^o,  and  the  first  two  terms  in  coefficient 

by, si  are  a  linear  combination  of  xo  and  Yo  (Eq.  (11)).  These  terms  match  the  terms  in  Eq.  (3)  for 
ye.  The  many  extra  terms  in  by,so  and  by,si  are  distortions  on  the  first  order  ROE  from  second 
order  gravitational  effects.  Center  point  Ze  equals  zero  in  the  first  order  framework  but  Zc  is  non¬ 
zero  in  the  second  order  framework.  Consequently,  all  three  terms  of  coefficient  bz,so  are  the 
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distortion  on  the  first  order  ROE  Ze.  A  summary  of  these  center  point  distortions  on  the  first 
order  ROE  is  given  in  Eq.  (12)  where  the  new  terms  are  not  explicitly  written  out  but  have 
obvious  definitions  from  the  given  data. 


x(t) 

1 

y<0 

=  {4-3Cq^x^xo 

1 

=  {6{Sn^x-*io'')>*0  +  <l>yo 

1 
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+  {|^-10C„^x+3C2„^x+12ii„xS„^x-12n^2)}xg 

1 

+  {|]^40Sn^x+3S2„j^x-22n„T-24noTC„^x)>*0 

1 

+  {3^Sn^x-no'*)>yo 
"“0 
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"“0 
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1 
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“3^0 

1 

“3^0 

1 
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1 
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1 
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+2{ii^jj^(-344C„^x-C2n„T)}Zo2o 

1 
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(10) 


x(t)  =  XcCt)  +  XhlCt)  +  Xshl(x)  +Xh2(T) 

y(t)  =  yc(i:)  +  yhiW  +  yshiW  +yh2(T:) 

z(t)  =  Zc  +  Zh](T)  +  Zshl(x)  +Zh2(x) 

where 

Xc('')  =  bx4io(l)  +  l>xjil(“o'')  +fex32(“o'f)^ 

yc(T)  =  by^(l)  +  by^iCiioX) 

Zc  =  bz^(l) 

Xhi(x)  =  axj,]Cos{0x4ii(T)}  ,  Xshi(x)  =ax^l(x)cos{0x^l(x)}  ,  xi^(x)  =ax^2cos{0x4i2(x)} 

yhl(T:)  =  ayaiisin{0yj,i(x)}  ,  yshl('')  =ay3hl(^)sin{0y,shl(^)}  .  yb2(^)  =ayji2si“-{0yji2('')} 

ZhlW  =  ’  Zshl('')  =az^hl(^)sin{0z^hl(''»  >  Zh2(^)  =azji2sin{®z4i2('')} 

=  “o’'  ■  ‘I’xhl  ’  ®x^hl('')  =  “o’'  ■  <l>x^hl  »  ®x4i2(''^)  =  ^HoX  -  c|>x^2 

0yjil(i:)  =  HoX  +  <|)yjji  ,  0y^hl('^)  =  +  <t>y^hl  >  »y4i2('')  =  +‘l>yji2 

~  ^0"®  “  ^qX  +  <t^z^2  ’  ®zji2^''^)  ^HqX  ■t‘*|*zj]2 


axhl  =  xjil  +  +  Cx^iKnoX)^ 

2  2  2  2  2  2  2 
ayjil  “  ^  yjil  ■•■  ^yjll  »  ^y^jjjCx)  =  {by^jjj  +  Cy^jj^Kn^x) 

^zjil  ~^zhl  ■*"  *'zjil  ’  ~  ■f^z^hl  ■*■  ^z^hlK^oX) 


2  u2  2 
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2.22 
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The  distortions  can  be  interpreted  as  additional  position,  velocity,  and  acceleration  sources. 
These  distortions  improve  the  characterization  of  the  deputy  motion  physics  including  curved  in- 
the-plane  departure  tracks  and  unbalanced  out-of-plane  harmonics  [31,32],  Note  if  all  second 
order  initial  condition  terms  are  set  to  zero,  all  distortions  become  zero  and  the  second  order 
ROE  collapse  to  the  first  order  ROE. 

The  y  axis  amplitude  2ae  is  constant  in  first  order  theory  (Eq.  (3))  and  also  in  second  order 
theory  for  the  first  harmonic  term  (Eq.  (10)),  but  is  statically  altered.  The  linear  initial  condition 
terms  xo,^o,y()  appearing  in  coefficients  by, hi  and  Cy,hi  (Eq.  (11))  after  being  squared  precisely 
form  the  first  order  4ae^  expression.  Other  second  order  initial  condition  terms  in  by, hi  and  Cy,hi 
represent  distortion  on  the  first  order  ROE  2ae.  In  a  similar  way,  the  linear  terms  zo,2o  in 
coefficients  bz,hi  and  Cz,hi  (Eq.  (11))  reproduce  the  first  order  z  axis  amplitude  be  (Eq.  (3))  and  the 
remaining  initial  condition  terms  cause  a  static  distortion.  Now  consider  a  comparison  of  y  axis 
harmonic  angles.  Both  first  order  angle  a  and  second  order  angle  0y,hi  have  the  same  frequency 
term  nox,  so  the  only  difference  for  a  and  0y,hi  originates  with  phase  angles  ao  =  a(0)  and  (t)y,hi. 
The  linear  initial  condition  terms  xo,Xo,yo  appearing  in  coefficients  by,hi  and  Cy,hi  (Eq.  (11)) 
exactly  match  the  terms  of  the  arctangent  argument  ratio  appearing  in  ao.  All  other  second  order 
initial  condition  terms  in  by,hi  and  Cy,hi  represent  distortion  on  the  phase  angle  of  the  first  order 
ROE  a.  Similarly,  the  linear  terms  zo,Zo  in  coefficients  bz,hi  and  Cz,hi  (Eq.  (11))  reproduce  the 
first  order  z  axis  phase  angle  Po  =  P(0)  (Eq.  (3))  and  the  remaining  initial  condition  terms  cause  a 
static  distortion.  A  summary  of  these  amplitude  and  phase  distortions  on  the  first  order  ROE  is 
given  in  Eq.  (13)  where  the  new  terms  are  not  explicitly  written  out  but  have  obvious  definitions 
from  the  given  data. 
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The  distortions  can  be  interpreted  as  additional  amplitude  and  phase  sources.  These  distortions 
also  improve  the  characterization  of  the  deputy  motion  physics  including  travel  distances  and 
travel  sequencing  in  both  in-plane  and  out-plane  axes  [31,32].  Note  if  all  second  order  initial 
condition  terms  are  set  to  zero,  all  distortions  become  zero  and  the  second  order  ROE  collapse  to 
the  first  order  ROE. 

Other  comparisons  can  be  made  between  the  first  order  and  second  order  ROE,  and  their 
impact  on  the  motion  characteristics.  Under  first  order  theory,  the  in-plane  harmonic  motion 
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maintains  a  fixed  yxx  =  2x1  amplitude  behavior,  regardless  of  the  initial  eonditions.  To  analyze 
the  corresponding  behavior  for  second  order  theory,  compare  the  expressions  for  amplitudes  ax, hi 
and  ay, hi,  or  bx,hi,  Cx,hi  and  by,hi,  Cy,hi.  For  every  initial  condition  term  in  bx,hi,  there  is  a 
corresponding  term  of  twice  the  value  but  with  opposite  sign  appearing  in  Cy,hi  with  two 
exceptions:  1)  the  terms  have  the  same  sign  and  no  doubling,  and  2)  Cy,hi  contains  the  extra 
term  yo^o.  Further,  although  two  terms  (^o,  zo^o)  in  Cx,hi  also  appear  in  by,hi  with  a  factor  of  2 
and  consistent  sign,  other  terms  are  not  scale  and  sign  consistent.  In  general  under  second  order 
theory,  the  in-plane  motion  does  not  strictly  follow  a  2-to-l  relationship,  although  an 
approximate  2-to-l  behavior  is  likely  due  to  the  many  matching  terms.  An  interesting 
observation  noted  in  Ref  [31]  is  that  two  of  the  secular-harmonic  coefficients  have  equal 
magnitudes  and  signs,  or  bx,shi  =  Cy,shi,  while  another  pair  of  the  secular-harmonic  coefficients 
have  equal  magnitudes  but  opposite  signs,  or  by,shi  =  -Cx,shi. 

Distortions  on  the  first  order  ROE  noted  previously  originate  with  the  center  point 
polynomials  and  first  harmonic  ellipses,  and  can  be  interpreted  as  parametric  distortions  where 
Eqs.  (12)-(13)  explicitly  describe  the  distortion  structure  and  how  it  affects  the  relative  motion. 
Other  distortions  exist  and  originate  with  the  first  secular-harmonic  and  second  harmonic 
ellipses.  These  effects  can  be  interpreted  as  coordinate  distortions  since  they  are  kept  as  separate 
terms  in  the  expanded  form  of  second  order  ROE.  If  the  coordinate  distortions  are  combined  and 
labeled  as  5x,  5y,  5z,  then  the  overall  comparison  of  second  order  and  first  order  ROE 
expressions  would  be  of  the  form 

A(0  =  e)}  +  a*j,i{ae^(a?)}cosrttoT-4i^j,l{ao^(ao)yS<a’i)}l  +  6x{x^i(-c),Xh2(t)> 

y(0  +ay4ii{ae^(4a^}sinLiiaX-t+yj,i{ao,5(a'o)^(aJ^}J+5y{ysjii{T)j'ij2(T)>  (14) 

z(t)  =  ZcfO  j8(Za)}  +  a5jj,i  {b  e  Wb^}sin[iioX+tz^i  {P  o  WP  o)»6(P  o)}]  +  6z{zRhl 

The  overscript  symbol  appearing  in  Eq.  (14)  distinguishes  two  different  parameter  sets  based 
on  different  trigonometric  forms  used  for  x(t)  in  the  first  order  Eq.  (3)  and  second  order  Eq.  (10) 
expressions.  The  Eq.  (14)  expressions  could  be  useful  in  performing  error  analysis  on  guidance, 
navigation,  and  control  algorithms  that  employ  first  order  dynamics,  or  exposing  critical  error 
sources  when  propagating  the  natural  relative  motion  using  first  order  dynamics. 

Eq.  (10)  indicates  overall  deputy  relative  orbital  motion  consists  of  superpositioning  (linear 
combining)  three  in-the-plane  and  three  out-of-plane  ellipses  and  polynomial  center  point  curves. 
Each  elliptic  motion  component  is  governed  by  a  set  of  three-dimensional  parametric  Lissajous 
functions  and  underlying  geometry  characterized  by  the  corresponding  six  relative  orbital 
elements  and  associated  auxiliary  circles  (see  Section  4),  giving  3x6  =  18  ROE  for  the  elliptic 
motion.  Movement  of  the  ellipse  centers  are  characterized  by  three  additional  ROE,  yielding 
18+3  =21  ROE  in  totality  using  the  expanded  formulation.  In  this  formulation,  two  of  the  ellipse 
sets  correspond  to  fixed  Lissajous  curves  and  the  other  ellipse  set  corresponds  to  variable 
amplitude  Lissajous  curves.  Component  motions  can  be  thought  of  as  natural  modes  that 
combine  to  form  overall  motion.  This  perspective  was  considered  in  Ref  [31]. 

To  further  explore  component  motions  and  the  underlying  geometry  described  by  second 
order  ROE,  a  three-dimensional  numeric  example  is  offered.  Deputy  initial  position  is  along  the 
LVLH  X  axis  ("above"  the  chief)  with  initial  velocity  components  in  all  three  x,  y,  and  z 
directions.  Circular  chief  orbital  parameters  and  deputy  relative  initial  conditions  are  specified  in 
Eq.  (15),  and  Fig.  13  shows  the  corresponding  relative  trajectory,  a  corkscrew  motion 
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representative  of  elose  proximity  eneounters,  for  approximately  2.5  orbits  where  T  denotes  chief 
orbital  period.  The  deputy  satellite  is  drifting  significantly  away  from  and  behind  the  chief 
satellite  while  simultaneously  oscillating  about  the  reference  orbit  plane.  After  two  orbits,  deputy 
position  is  approximately  -800  km  along-track  and  -40  km  radial  with  ±20  km  cross-track 
amplitudes  during  the  motion.  Two  trajectories  are  shown  in  Fig.  13:  the  "exact"  trajectory 
generated  by  nonlinear  simulation  (NTS)  based  on  Runge-Kutta  4*  order  numerical  integration 
and  the  approximate  QV  trajectory  from  Eq.  (9).  The  overlay  plot  implies  orbit  propagation  error 
due  to  differences  between  two-body  exact  and  approximate  QV  motion  models  will  be  minimal 
for  the  temporal  horizon  shown.  Simulation  parameters  include  Earth  gravitation  constant  p  = 
398,600  km^/s^  and  integration  time  step  At  =  (27r/no)/3600  with  rad/s  units  for  no. 

Chief:  Deputy:  to  =  Os 

Rj,  =  7100  km  Xf,  =  0.2  km  ,  yQ  =  0  km  ,  =  0  km  (15) 

Hq  =  0.001055  rad/s  Xq  =  0.002  km/s  ,  yg  =  0.02  km/s  ,  Zg  =  0.02  km/s 

The  relative  trajectory  for  8  orbits  (8T)  is  shown  in  Fig.  14  and  is  used  to  investigate  the 
motion  geometry  and  second  order  ROE  based  on  the  expanded  formulation.  Four  specific  times 
ti  =  1/8T,  t2  =  1/4T,  t3  =  IT,  U  =  2T  are  highlighted  along  the  trajectory  and  used  in  the 
discussion.  Fig.  15  shows  the  four  component  trajectories:  secular  center  point  polynomial 
motion  (c),  first  harmonic  elliptic  motion  (hi),  first  secular-harmonic  elliptic  motion  (shl),  and 
second  harmonic  elliptic  motion  (h2).  Each  curve  shows  the  start  and  stop  positions,  and  Tab.  1 
lists  the  corresponding  initial  and  final  values  for  the  twenty-one  ROE.  Note  each  three- 
dimensional  elliptic  curve  can  be  projected  to  the  xy  and  xz  planes  (shown  in  Fig.  15)  with 
corresponding  auxiliary  circles  and  deputy  shadow  satellites  (not  shown  in  Fig.  15)  to  visualize 
and  quantify  the  motion  geometry.  For  example,  the  four  start  points  in  Fig.  15  can  be 
graphically  superimposed  to  form  the  overall  start  point  in  Fig.  14,  or  mathematically  for  x(to) 


Table  1.  ROE  Initial  and  Final  Values  -  Expanded  Formulation 


Units 

km 

km 

km 

deg 

deg 

deg 

Initial  Val. 

ROE 

Xc 

yc 

Zc 

Center  Pt. 

3.8899e+01 

-3.8059e+00 

7.5881e-03 

ROE 

ax,hl ,  ax,  shl  ,ax,h2 

ay  ,h  1 ,  ay ,  sh  1 ,  ay  ,h2 

az  ,h  1 ,  az ,  sh  1 ,  az  ,h2 

0x,hl  ,0x,shl  ,0x,h2 

0y,hl  ,0y,shl  ,0y,h2 

0z,hl,0z,shl,0z,h2 

1  St  Har. 

3.8864e+01 

7.7727e+01 

1.9004e+01 

-1.7715e+02 

2.8028e+00 

-3.0503e-02 

1st  See. -Har. 

0 

0 

0 

-8.7182e+01 

-8.7182e+01 

-9.0000e+01 

2nd  Har. 

1.1725e-01 

3.9762e-02 

5.1450e-02 

5.0288e+00 

7.4257e+00 

2.8179e+00 

Final  Val. 

ROE 

Xc 

yc 

Zc 

Center  Pt. 

-5.6236e+02 

-2.9417e+03 

7.5881e-03 

ROE 

ax,hl ,  ax,  shl  ,ax,h2 

ay  ,h  1 ,  ay ,  sh  1 ,  ay  ,h2 

az  ,h  1 ,  az ,  sh  1 ,  az  ,h2 

0x,hl  ,0x,shl  ,0x,h2 

0y,hl  ,0y,shl  ,0y,h2 

0z,hl,0z,shl,0z,h2 

1  St  Har. 

3.8864e+01 

7.7727e+01 

1.9004e+01 

2.7028e+03 

2.8828e+03 

2.8800e+03 

1st  See. -Har. 

1.5844e+01 

1.5844e+01 

7.7894e+00 

2.7928e+03 

2.7928e+03 

2.7900e+03 

2nd  Har. 

1.1725e-01 

3.9762e-02 

5.1450e-02 

5.7650e+03 

5.7674e+03 

5.7628e+03 

Approved  for  public  release;  distribution  is  unlimited. 

24 


2  tkm)  2  iV.m) 


Figure  13.  QV  vs.  NLS  for  2.5T  Figure  14.  QV  for  8T  -  Expanded  Formulation 
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Figure  15.  Component  Motions  -  Expanded  Formulation 
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x(to)  =  X^(to)  +  Xhjfto)  +  XgijiCto)  +  Xj^(to) 

Xo  =bx^  +  ajijiiCO»{e*j,i(to)}  +  ax^bl(to)«»<0x^l(W}  +  axJi2CO8{0x^(to)} 

0.2  km  =38.899  +  38.864cos{-177.15*}  +  0cos{-87.182"}  +0.11725cos{5.0288“}km  (16) 

0.2  km  =38.899  + -38.8159  +  0  +  0.11680  km 

0.2  km  =0.1999  km  V  (0.1  mTiuncatioB  Einjr) 

Although  a  similar  process  can  be  carried  out  for  y(to),  z(to),  or  more  generally  for  x(t),  y(t),  z(t) 
at  any  specific  time  t,  emphasis  is  now  directed  to  higher  level  geometric  characteristics. 

Several  observations  are  made  from  Fig.  15.  The  center  point  curve  consists  of  quadratic 
propagation  along  the  x  axis  and  linear  propagation  along  the  y  axis,  while  the  whole  curve  is 
offset  above  the  z  =  0  km  plane  by  a  small  constant  amount  (zc  =  7.6  m,  see  Tab.  1),  causing  the 
center  point  curve  to  lie  in  a  plane  parallel  to  the  chief  orbital  plane.  The  first  harmonic  curve 
consists  of  a  closed  fixed  Lissajous  ellipse  lying  in  a  plane  skewed  to  the  chief  orbital  plane.  The 
first  secular-harmonic  curve  consists  of  an  open  variable  amplitude  Lissajous  elliptical  spiral, 
also  lying  in  a  skewed  plane.  Finally,  the  second  harmonic  curve  also  corresponds  to  a  closed 
fixed  Lissajous  ellipse  lying  in  a  skewed  plane.  Note  in  Fig.  15  the  center  point  x,y  amplitudes 
are  at  least  an  order  of  magnitude  larger  than  for  the  first  harmonic  amplitudes,  in  turn,  the  first 
harmonic  x,y,z  amplitudes  are  approximately  one  magnitude  order  larger  than  the  first  secular- 
harmonic  amplitudes,  and  in  turn  again,  the  second  harmonic  amplitudes  are  much  smaller  than 
the  first  secular-harmonic  amplitudes.  In  other  words,  the  in-plane  center  point  components 
determine  to  a  large  extent  the  overall  nature  of  the  deputy  motion,  and  the  first  harmonic 
component  introduces  a  three-dimensional  perturbing  oscillation  about  the  parabolic  path. 
Likewise,  the  first  secular-harmonic  components  behave  like  a  three-dimensional  expanding 
oscillation  perturbation  about  the  first  harmonic  path,  and  all  of  the  above  is  perturbed  by  a 
higher  frequency  oscillation  originating  from  the  second  harmonic  component. 

When  viewed  from  above  the  z  =  0  km  plane,  the  first  harmonic  motion  is  in  an  opposite 
direction  to  the  first  secular-harmonic  and  second  harmonic  motions  (clockwise  for  hi,  counter¬ 
clockwise  for  shl,  h2).  Careful  examination  of  the  sequencing  of  indexed  points  ti,t2  for  the  hi 
and  h2  ellipses  in  Fig.  15  illustrates  the  opposing  motions.  Such  peculiar  behavior  is  a  trait 
exhibited  by  Lissajous  curves  under  certain  parameter  conditions.  When  the  x,y  or  x,z  angular 
ROE  become  sufficiently  out-of-phase,  the  direction  of  motion  along  the  Lissajous  path  can 
reverse  (0x,hi  =  -177  deg  vs.  0y,hi  =  2.8  deg  or  0z,hi  =  -0.031  deg,  see  Tab.  1).  The  variable  phase 
Lissajous  curve  example  in  Fig.  5  for  the  0  <  x  <  4(27r)  case  exhibited  similar  behavior  when  the 
time  dependent  phase  reached  a  critical  value.  Also,  for  every  increment  of  first  harmonic  (and 
first  secular-harmonic)  angular  motion  occurring  over  a  unit  of  time,  approximately  "two" 
increments  of  second  harmonic  angular  motion  simultaneously  occur.  Careful  examination  of  the 
spacing  of  indexed  points  ti,t2  for  the  hi  and  h2  ellipses  in  Fig.  15  illustrates  the  differing 
increments.  This  behavior  is  a  consequence  of  the  nox  vs.  2noX  frequencies.  Opposing  travel 
directions  and  differing  travel  rates  in  the  various  elliptic  motion  components  leads  to  quadrant 
dependent  perturbations  that  push  the  deputy  satellite  outside  of,  or  pull  the  deputy  satellite 
inside  of,  the  fundamental  harmonic  elliptical  path.  Fig.  16  shows  the  elliptical  component 
motions  overlaid.  Note  the  required  shl  and  h2  scalings  (x2.5  for  shl  at  t  =  8T,  x250  for  h2) 
required  to  make  the  corresponding  curves  equally  visible  to  the  hi  curve.  The  side-view 
perspective  shows  that  each  plane  containing  the  Lissajous  ellipses  have  different  relative 
"ascending  longitudes"  and  "inclinations",  further  enriching  the  geometry  underlying  the  overall 
deputy  motion. 
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Figure  16.  Overlaid  Elliptic  Component  Motions  -  Expanded  Formulation 

Fig.  17  shows  the  twenty-one  second  order  ROE  time  evolutions  for  the  example.  Note  the 
radial  and  cross-track  center  point  ROE  are  scaled  (x2.5  for  Xc,  x25000  for  Zc)  to  make  the 
corresponding  response  curves  equally  visible  to  the  along-track  curve  (yc).  The  quadratic  and 
linear  drifts  and  constant  offset  behaviors  in  the  center  point  coordinates  are  clearly  noticeable. 
Observe  how  the  first  harmonic  and  second  harmonic  ROE  amplitudes  are  constant  (see  Tab.  1 
for  numeric  correlation),  but  the  first  secular-harmonic  ROE  amplitudes  are  linearly  increasing 
from  the  initial  values  (see  Tab.  1).  For  the  hi  component,  note  the  amplitude  ratio  ay,hi/ax,hi  = 
1.99997  is  not  precisely  equal  to  2  but  is  practically  2  as  discussed  in  the  analysis  of  coefficients 
bx,hi,  Cx,hi  and  by, hi,  Cy,hi.  For  the  shl  component,  the  ratio  of  amplitude  rate  of  change  is 
d(ay,shi)/dt/d(ax,shi)/dt  =  1,  and  for  the  h2  component,  the  amplitude  ratio  is  ay,h2/ax,h2  =  0.339.  The 
shl  amplitude  rate  ratio  equaling  1  is  a  consequence  of  the  noted  conditions  by, shl  ~  -Cx,shl  and 
bx,shi  =  Cy,shi  discussed  previously.  In  all  elliptical  motions,  the  ROE  angles  are  linearly 
increasing  from  their  initial  values  (see  Tab.  1).  In  Fig.  17,  the  x,y  ROE  angles  for  the  shl 
component  are  precisely  equal  for  all  times  or  initially  0x,shi  =  -87.182  deg  =  0y,shi  (see  Tab.  1), 

again  as  a  consequence  of  conditions  by, shl Cx,shl  and  bx,shi  =  Cy,shi.  Fig.  17  gives  the  impression 

that  0z,shi  for  shl  is  also  equal  to  this  value,  however,  0z,shi  is  slightly  different  due  to  the  initial 
phase  offset  (0z,shi  =  -90.000  deg,  see  Tab.  1).  For  the  h2  component,  all  three  ROE  angles 
appear  equal  but  differ  slightly  due  to  the  small  dissimilar  initial  phases  (see  Tab.  1).  These  small 
phase  differences  for  shl  and  h2  correlate  to  the  counter-clockwise  or  direct  movement  around 
the  respective  ellipses.  In  contrast,  the  hi  angular  ROE  Ox, hi  evolution  in  Fig.  17  is 
approximately  180  deg  behind  the  0y,hi  and  0z,hi  responses,  causing  the  clockwise  or  retrograde 
movement  around  the  first  harmonic  ellipse  noted  previously.  This  example  has  demonstrated 
how  the  expanded  formulation  of  second  order  relative  orbital  elements  characterize  the 
underlying  geometry  of  the  deputy  satellite  relative  motion.  This  characterization  may  be  useful 
for  formulating  strategy  for  predictive  propagation,  maneuvering  planning,  rendezvous-intercept, 
or  determination-estimation. 
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Figure  17.  Second  Order  Relative  Orbital  Element  Responses  -  Expanded  Formulation 
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4.2  SECOND  ORDER  RELATIVE  ORBITAL  ELEMENTS  -  COMPACTED  FORM 


An  alternate  formulation  of  second  order  ROE  is  possible  by  expressing  the  nonlinear  closed- 
form  QV  relative  trajectory  expressions  for  x(t),  y(t),  and  z(t)  in  a  different  manner.  By  using 
trigonometric  identities  to  replace  second  harmonics  with  powers  of  first  harmonics,  grouping  of 
terms  in  Eq.  (9)  can  be  accomplished  on  the  individual  expansion  functions:  1,  nox,  (nox)^, 
cos(nox),  sin(nox),  five  functions  rather  than  nine  in  the  expanded  formulation.  Expressions  in  Eq. 
(9)  can  be  converted  to  the  form  in  Eq.  (17),  where  Eqs.  (18)  and  (11)  provide  lower  level 
coefficient  definitions.  Eq.  (17)  indicates  trajectory  expressions  contain  one  set  of  three- 
dimensional  parametric  Lissajous  functions  for  the  effective  first  harmonic  noX  indicated  by 
subscript  "hie",  in  addition  to  a  set  of  polynomial  functions  of  noX  indicated  by  subscript  "c". 
The  Lissajous  functions  produce  one  elliptical  space  curve  and  correspondingly  one  set  of  xy  and 
xz  planar  ellipses,  and  the  polynomial  functions  produce  a  secular  quadratic  space  curve. 
Following  Sections  3-4,  a  set  of  nine  ROE  are  defined  directly  from  the  compacted  solution  form 
in  Eq.  (17)  and  include:  {xc,  yc,  Zc}  and  {ax, hie,  ay, hie,  az,hie,  0x,hie,  0y,hie,  0z,hie}.  Elements  (xc,  yc, 
Zc}  represent  the  in-plane  and  out-plane  center  point  coordinates  for  the  two-dimensional 
ellipses.  Elements  (ax, hie,  ay, hie,  az,hie,  0x,hie,  0y,hie,  0z,hie}  characterize  the  size-shape  of,  and 
angular  position  along,  the  two-dimensional  ellipses  corresponding  to  effective  harmonic  noX. 
Recall  each  ellipse  will  have  two  associated  auxiliary  circles  that  are  explicitly  described  by  the 
elements,  which  in  turn  implicitly  describe  the  actual  geometry  of  the  motion  along  the  ellipse  in 
a  meaningful  way.  In  this  formulation,  all  three  amplitudes  and  all  three  phases  are  time 
dependent  and  correspond  to  a  variable  Lissajous  curve.  The  two  different  formulations, 
expanded  formulation  in  the  previous  section  and  the  compacted  formulation  in  this  section, 
provide  a  clear  trade-off  regarding  complexity:  number  of  relative  orbital  elements  vs.  time 
dependent  amplitudes-phases. 


x(t)  =  Xc(t)  +Xiiie(T) 

y(t)  =yc(t)+yhie('*)  (17) 

z(t)  =  Zg  +Zhle('') 
where 

Xc(T)  =bx^l)  +  bx;»l(llo‘<^)+b;i^{not)^  , 

yc(T)  =  by_jo(l)  +  by^i(noT)  ,  yhle(^)  =«yjile('')®“<0y4ile('')> 

z®  =  1>Z40C1)  •  =azJile('t)«Q{0zJileW} 

"t"  tyjile^^^  »  ayjilo^^)  “  ^yjile(”®^  »  ^aa'{0y^jg(T)}' =  byjjjg(T) 

=  “o'*  +  tzJileC*)  »  “zJilcW  =  I’zJileW  +  t^zJileW  »  =  l*z4ile(*)  /^ziileC*) 
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^  **■  ^  x,shl(®o'^)  ■*"  ^xJjIc^DqT  ■*■  ^iJila^npT 

■*“  '*’  ^aJiIc^iIqT  "*■  ^'xJils^npT 

^y,hle(^)  "*■  '*’ ^y^ilc^n^T: ■*‘^y4ils^iipX 

*'y^ile(‘^)  ~*'y4il(l)  +  *'y^l(**o"'')  ■*"  ^yJilc^ttQX  ■*"  ^y^ils^Qo'* 
^  z4ile(^)  ■*■  ^  z^lC^oi^)  "*"  ^zJilc^DoX  '*' ^ zJiIs^iIqT 

^zjiieW  ■*■  "*■  *'zJilc^nQX  ■*"  ^zJiIs^DqX 
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L  g  *0*0  ■  2k  *0^0 

^DoRo  njEo 

-ZQyo+2  2^  ^ozq 
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‘o^o  njRp 

^Z0*0  +  -^^0i0 


(18) 


A  comparison  of  second  order  and  first  order  ROE  expressions  is  considered  next. 
Expressions  for  center  point  coordinates  Xc,  yc,  Zc  are  identical  for  both  the  compacted  and 
expanded  formulations  (see  Eqs.  (17)  vs.  (10)),  where  coordinate  Xc  is  time  dependent  and 
involves  secular  terms  of  powers  0,1,2,  coordinate  yc  is  time  dependent  and  involves  secular 
terms  of  powers  0, 1 ,  and  coordinate  Zc  is  a  constant.  Further,  the  lower  level  secular  polynomial 
coefficients  with  subscripts  "sO",  "si",  "s2"  remain  unchanged  from  the  Eq.  (11)  definitions  and 
are  tunctions  of  the  initial  positions  and  velocities.  Therefore,  Eq.  (12)  describing  distortions  on 
the  first  order  center  point  coordinate  ROE  also  applies  here.  With  the  secular-harmonic  terms 
and  higher  frequency  terms  absorbed  into  the  fundamental  harmonic  terms,  the  compacted 
formulation  provides  a  cleaner  comparison  between  the  first  order  and  second  order  amplitude 
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and  angle  ROE.  The  y  axis  amplitude  2ae  is  eonstant  in  first  order  theory  (Eq.  (3))  but  the 
eorresponding  amplitude  ay,hie(T)  for  the  first  harmonic  effective  term  is  time  dependent  (Eq. 
(17)).  Einear  initial  condition  terms  xo,^(),yo  appearing  in  coefficients  by,hi  and  Cy,hi  within 

by, hie(T)  and  Cy,hie(T)  (Eqs.  (1 1)  and  (18))  form  the  first  order  4ae^  expression  after  being  squared. 
All  other  second  order  initial  condition  terms  within  by,hie(T)  and  Cy,hie(T)  originating  from 
coefficients  by, hi,  by.shi,  by,hic,  by, his  and  Cy,hi,  Cy,shi,  Cy,hic,  Cy,his  augment  the  y  axis  amplitude 
statically  and  dynamically.  The  dynamic  distortion  involves  a  linearly  increasing  term  and 
harmonic  oscillating  terms.  In  a  similar  way,  the  linear  terms  20,^0  in  coefficients  bz,hi  and  Cz,hi 
within  bz,hie(T)  and  Cz,hie(T)  (Eqs.  (11)  and  (18))  generate  the  first  order  z  axis  amplitude  be  (Eq. 
(3)).  All  remaining  second  order  initial  condition  terms  within  bz,hie(T)  and  Cz,hie(T)  coming  from 

bz, hi,  bz,shi,  bz,hic,  bz,his  and  Cz,hi,  Cz,shi,  Cz,hic,  Cz.his  modulatc  the  z  axis  amplitude.  The  modulation 
again  involves  a  static  effect,  a  secular  effect,  and  a  harmonic  effect.  Therefore,  relations 
between  first  order  ae,  be  and  second  order  ay,hie(T),  az,hie(T)  are  similar  to  Eq.  (13)  but  with 
different  implied  function  definitions  for  the  distortions.  Moving  on  to  the  angular  coordinate 
comparison,  the  only  difference  between  the  first  order  and  second  order  ROE  lies  in  the  phase 
angles.  For  the  y  axis,  the  comparison  is  between  first  order  phase  angle  ao  and  second  order 
phase  angle  <|)y,hie(T).  Linear  initial  condition  terms  xo,^o,yo  that  generate  ao  (see  Eq.  (3))  can  be 
found  within  coefficients  by,hie(T)  and  Cy,hie(T)  that  contribute  to  (|)y,hie(T)  (Eqs.  (17)  and  (18)).  All 
other  second  order  initial  condition  terms  within  (t)y,hie(T)  are  static  offset  or  time  dependent 
distortions  to  ao.  Similar  conclusions  can  be  drawn  for  Po,  (|)z,hie(T)  and  zo,^o.  Relationships 
between  first  order  ao,  Po  and  second  order  (|)y,hie(T),  (|)z,hie(T)  are  similar  to  Eq.  (13)  but  with 
different  implied  function  definitions  for  the  distortions.  In  summary,  Eq.  (19)  describes  the 
parametric  distortions  for  the  compacted  formulation  of  second  order  ROE.  Note  the  absence  of 
direct  coordinate  distortions  (like  in  Eq.  (14))  since  all  extra  terms  have  been  folded  into  the 
effective  first  harmonic  term.  As  before,  if  all  second  order  initial  condition  terms  are  set  to  zero, 
then  all  distortions  become  zero  in  Eq.  (19),  and  the  second  order  ROE  reduce  to  the  first  order 
ROE.  Eq.  (19)  again  could  be  useful  in  specific  applications  concerning  the  control  and 
dynamics  of  relative  satellite  motion. 

Several  interesting  observations  regarding  effective  first  harmonic  coefficient  relationships 
are  made  here.  First,  previously  noted  relations  (Section  5)  between  the  in-plane  x  and  y  axis 
coefficients  are  also  true  here:  2bx,hi  -Cy,hi,  by,hi  2cx,hi  and  bx,shi  =  Cy.shi,  by,shi  =  -Cx,shi.  Second, 
many  additional  harmonic  terms  within  bx,hie(T)  and  Cx,hie(T),  by,hie(T)  and  Cy,hie(T),  and  within 
bz,hie(T)  and  Cz,hie(T),  are  equivalent.  These  additional  relations  are  summarized  as  bx,hic  =  -Cx,his, 
bx.his  “  Cx.hic,  by,hic  “  -Cy.his,  by.his  “  Cy.hic,  bz.hic  “  -Cz.his,  bz.his  “  Cz,hic.  Only  approximate 
conclusions  can  be  drawn  from  these  relations.  The  y-to-x  amplitude  behavior  will  initiate  with 
an  approximate  2-to-l  relation,  the  secular  amplitude  growth  between  the  y  and  x  axes  will 
follow  an  approximate  1-to-l  relation,  and  although  the  harmonic  oscillation  amplitudes 
associated  with  a  given  axis  are  balanced  between  the  cosine  and  sine  terms,  the  y-to-x  harmonic 
amplitude  relation  does  not  necessarily  follow  an  approximate  2-to-l  or  1-to-l  rule.  Eq.  (17) 
indicates  the  overall  deputy  relative  orbital  motion  consists  of  the  direct  sum  of  the  in-plane  and 
out-plane  ellipses  and  polynomial  center  point  curves.  The  elliptical  motion  is  characterized  by 
the  six  relative  orbital  elements  and  the  associated  auxiliary  circles  from  the  variable  Lissajous 
functions,  and  movement  of  the  ellipse  centers  are  characterized  by  three  additional  ROE, 
yielding  6+3  =  9  ROE  in  totality  using  the  compacted  formulation.  These  component  motions 
can  also  be  interpreted  as  natural  modes  that  combine  and  form  the  overall  motion. 
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A(t)  =  Xc{'t;teAXe)^(^e)^^e)>  +a;j4,ie{ae^(a5}cos[llo't-<|>3t^le{aoWo)Wo)}] 

y(t)  =  yc{'«,ye»^e*S(ye).S(:^e)}  +ayj,ie{ae^4ae)}sm[llo-t4<tiyjiie{ao^(ai,)^(ao)}]  (19) 

Z(t)  =  Zg{O^Ze)}  +a2^1e{be^^}sm[noT+<|)2^1e{Po»6(Po)APo)}l 


fyjile  =  »  tan(<l>yjiie)  =  [g>n(«o)  +  taP{&(ao)}I  /  [1  +  ^{&(ao)}] 

2Dd  Order  Ist  Order  DlstortioDS  2iid  Order  let  Order  Dishxtioiis  Distoidons 

az^ile  =  ^  +  5^^  »  tan(<|>2jjie)=[tm^  +  tan{aOo)}]/[l+taii{&(po)}I 

2Dd  Order  1*^  Order  DistortioDS  '"^dOidtt^  1*^  Order  DistoitioiiB  Distoitioiis 

The  previous  numeric  example  in  Eq.  (15)  is  revisited  here  to  further  explore  the  component 
motions  and  underlying  geometry  described  by  the  alternate  second  order  ROE.  The  relative 
trajectory  for  0  <  t  <  8T  is  shown  in  Fig.  18  and  is  used  to  examine  the  motion  geometry  and 
second  order  ROE  based  on  the  compacted  formulation.  Two  specific  times  ts  =  (l+3/8)T,  te  = 
(3+3/8)T  are  highlighted  along  the  trajectory  and  used  in  the  discussion.  Fig.  19  shows  the  two 
component  trajectories:  secular  center  point  polynomial  motion  (c),  and  effective  first  harmonic 
elliptic  motion  (hie).  Each  curve  shows  the  start  and  stop  positions,  and  Tab.  2  lists  the 
corresponding  initial  and  final  values  for  the  nine  ROE.  Note  the  three-dimensional  elliptic  curve 
can  be  projected  to  the  xy  and  xz  planes  (shown  in  Fig.  19)  with  corresponding  auxiliary  circles 
and  deputy  shadow  satellites  (not  shown  in  Fig.  19)  to  visualize  and  quantify  the  motion 
geometry.  For  example,  the  two  start  points  in  Fig.  19  can  be  graphically  superimposed  to  form 
the  overall  start  point  in  Fig.  18,  or  mathematically  for  y(to) 

y(to)  =  yc<to)  +  yhie(to) 

yo  =  byjrf)  +  ay_hie(to)sm{0yj,je(to)} 

0  km  =  -3.8059  +  77.767am{2.8052‘*}  km  (20) 

0  km  =  -3.8059  +  3.8059  km 

0km=0km  V 

Although  a  similar  process  can  be  carried  out  for  x(to),  z(to),  or  more  generally  for  x(t),  y(t),  z(t) 
at  any  specific  time  t,  emphasis  is  now  directed  to  higher  level  geometric  characteristics. 
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Table  2.  R( 

3E  Initial  and 

Final  Values  -  Compactec 

Formulation 

Units 

km 

km 

km 

deg 

deg 

deg 

Initial  Val. 

ROE 

Xc 

yc 

Zc 

Center  Pt. 

3.8899e+01 

-3.8059e+00 

7.5881e-03 

ROE 

ax, hie 

ay,hle 

az,hle 

0x,hle 

0y,hle 

0z,hle 

1st  Har.  Eff. 

3.8747e+01 

7.7767e+01 

1.9055e+01 

-1.7716e+02 

2.8052e+00 

-2.2816e-02 

Final  Val. 

ROE 

Xc 

yc 

Zc 

Center  Pt. 

-5.6236e+02 

-2.9417e+03 

7.5881e-03 

ROE 

ax, hie 

ay,hle 

az,hle 

0x,hle 

0y,hle 

0z,hle 

1st  Har.  Eff. 

4.1867e+01 

7.9368e+01 

2.0589e+01 

2.7251e+03 

2.8713e+03 

2.8577e+03 

Figure  18.  QV  for  8T  -  Compacted  Formulation 
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Figure  19.  Component  Motions  -  Compacted  Formulation 
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Several  observations  are  made  from  Fig.  19.  The  eenter  point  curve  again  lies  in  a  plane 
slightly  above  the  chief  orbital  plane  (zc  =  7.6  m,  see  Tab.  2)  and  consists  of  quadratic 
propagation  along  the  x  axis  and  linear  propagation  along  the  y  axis.  The  first  harmonic  effective 
curve  consists  of  an  open  variable  Lissajous  ellipse.  This  space  curve  is  three-dimensional  and 
does  not  lie  in  any  fixed  plane.  Note  in  Fig.  19  the  center  point  x,y  amplitudes  are  at  least  an 
order  of  magnitude  larger  than  for  the  first  harmonic  effective  amplitudes.  In  other  words,  the  in¬ 
plane  center  point  components  determine  to  a  large  extent  the  overall  nature  of  the  deputy 
motion,  and  the  first  harmonic  effective  component  introduces  a  three-dimensional  perturbing 
oscillation  about  the  parabolic  path.  The  direction  of  first  harmonic  effective  motion  opposes  the 
counter-clockwise  buildup  of  angles  0x,hie  and  0y,hie,  or  angles  0x,hie  and  0z,hie,  in  the  xy  and  xz 
planes.  Careful  examination  of  the  sequencing  of  indexed  points  ts  and  te  from  the  start  point  in 
Fig.  19  illustrates  the  retrograde  motion.  The  large  phase  angle  difference  noted  in  Tab.  2 
between  the  x  and  y,z  axes  causes  the  retrograde  behavior  (0x,hie  =  -177  deg  vs.  0y,hie  =  2.8  deg 
and  0z,hie  =  -0.023  deg).  Note  the  similarity  of  the  first  harmonic  effective  component  motion  in 
Fig.  19  with  the  variable  amplitude-phase  Lissajous  space  curve  example  in  Section  4  for  the  0a 
=  +30  and  -100  deg  azimuth  views  in  Fig.  12.  In  particular,  the  three-dimensional  curve  of  Fig. 
19  exhibits  signatures  of  the  conchiglie  pasta  or  seashell  shape  and  multi-coil  spring  shape.  Also, 
the  xz  projected  curve  in  Fig.  19  is  highly  similar  to  the  example  Lissajous  planar  curve  with 
variable  amplitude  and  phase  in  Fig.  7  for  0  <  x  <  8(27r).  Although  not  visible  in  Fig.  19,  the 
elliptical  curve  also  experiences  small  oscillations  along  the  entire  length. 

Fig.  20  shows  the  nine  second  order  ROE  time  evolutions  for  the  example.  Note  the  radial 
and  cross-track  center  point  ROE  are  again  scaled  (x2.5  for  Xc,  x25000  for  Zc)  to  make  the 
corresponding  response  curves  equally  visible  to  the  along-track  curve  (yc).  The  quadratic  and 
linear  drifts  and  constant  offset  behaviors  in  the  center  point  coordinates  are  clearly  noticeable. 
The  time  dependent  nature  of  the  elliptical  motion  amplitudes  is  seen  in  Fig.  20,  although  over  8 
orbits  the  amplitudes  only  change  by  several  kilometers.  With  a  tighter  axis  view,  the  secular 
growth  and  harmonic  oscillation  in  the  ROE  amplitudes  is  visible  where  the  y  axis  oscillation 
amplitude  is  much  smaller  than  the  x  and  z  axis  behavior.  Further,  the  y  axis  secular  growth  lags 
the  corresponding  behavior  in  the  x  and  z  axes  by  a  small  amount.  The  y-to-x  amplitude  ratio  at 
the  initial  time  is  ay,hie/ax,hie  =  2.007  and  after  8  orbits  the  ratio  is  ay,hie/ax,hie  =  1.896.  Although 
the  precise  2-to-l  amplitude  nature  is  not  present  in  the  second  order  motion,  the  actual 
amplitude  relation  is  not  far  from  the  first  order  behavior,  as  predicted  during  the  examination  of 
coefficients  bx,hie(x),  Cx,hie(T)  and  by,hie(T),  Cy,hie(T).  The  z-to-x  amplitude  ratios  at  the  initial  and 
final  times  are  az,hie/ax,hie  =  0.492  and  0.492  even  though  the  ratio  varies  by  a  small  amount  over 
time.  In  the  elliptical  motion,  the  ROE  angles  are  linearly  increasing  from  their  initial  values  (see 
Tab.  2).  Fig.  20  gives  the  impression  that  angles  Oy.hie  and  0z,hie  are  equal,  however  these  angles 
are  slightly  different  due  to  the  initial  phase  offset  (0y,hie  =  2.8  deg  vs.  0z,hie  =  -0.023  deg,  see 
Tab.  2).  In  contrast,  the  angular  ROE  Ox, hie  evolution  in  Fig.  20  is  again  approximately  180  deg 
behind  the  0y,hie  and  0z,hie  responses,  causing  the  clockwise  or  retrograde  movement  around  the 
elliptical  space  curve  noted  previously.  By  plotting  only  the  phase  angles  and  removing  the  out- 
of-phase  trait  in  the  x  axis,  the  time  dependent  nature  of  the  phase  angles  in  each  axis  of  elliptical 
motion  can  be  observed  (see  Fig.  20).  This  example  has  demonstrated  how  the  compacted 
formulation  of  second  order  relative  orbital  elements  characterize  the  underlying  geometry  of  the 
deputy  satellite  relative  motion. 
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4.3  LISSAJOUS  CURVES  -  REVISITED 


Before  considering  second  order  conditions  for  periodic  relative  motion,  conditions  for  fixed 
and  variable  closed  Lissajous  curves  are  examined.  These  topics  are  closely  related  to  periodicity 
based  on  second  order  ROE.  Note  the  following  analysis  emphasizes  whole  period  return  points 
and  does  not  address  rational  sub-period  or  super-period  self-intersection  points.  First,  consider 
the  case  of  variable  frequency  but  constant  amplitude-phase,  or 

x(t)  =axCOs{w^(T)T-<t)x} 

y(t)  =aySin{a>y(T)T-Kl)y}  ^  ^ 

Parametric  conditions  are  sought  for  the  situation  where  coordinate  values  x,y  simultaneously 
return  to  previous  values  after  repeat  time  Tr  indicated  in  Eq.  (22).  Note  Tr  represents  a  change  in 
time,  i.e.,  Tr  =  At  or  Ax. 


x(t+Tr)  =x(t) 
y(t+T,)  =y(t) 

The  equivalent  angular  conditions  from  Eq.  (22)  after  canceling  out  equal  phase  terms  are 

oix('^+Tr)(t+Tr)  =(2n;n)nx 

^ _ _ ^ 

nxfx(T,Tr) 

0)y(T-|-Tr)(T-|-Tr)  -03y(T)T  =(2jin)ny 

' _ _ _ ^ 

nyfy(T,T^) 


(22) 


(23) 


where  n  represents  a  common  integer  number  of  whole  function  revolutions  over  time  Tr  while 
nx  and  ny  denote  revolution  proportionality  integers  between  the  x  and  y  axes.  For  a  specific  x, 
Eq.  (23)  represents  two  conditions  that  the  one  variable  Tr  must  simultaneously  satisfy  for 
repeating  x,y.  The  closed  curve  requirement  means  the  ratio  of  the  x  and  y  axis  frequency 
function  differences  evaluated  at  two  specific  times  separated  by  the  repeat  time  must  be  rational, 

Hx  _  m^(t-hT,.)(t-hT,.)  -  04) 

II  y  tOy(T-|-Tj.)(T-|-T,.)  -  tOy(T)T 

otherwise  the  curve  is  open.  With  general  time  x,  Eqs.  (23)-(24)  greatly  restrict  the  class  of  time 
dependent  frequency  functions  that  are  capable  of  producing  a  closed  curve.  One  class  of 
sufficient  but  not  necessary  frequency  functions  for  repeating  x,y  are  rationally  proportional 
cox(x)  and  coy(x),  or  cox(x)  and  coy(x)  are  obtained  from  whole  number  scaling  of  a  base  function 
co(x)  as  indicated  below. 


to^fx)  =nj(^a)(x) 

tt»y(x)  =nyO)(x) 

In  this  case,  the  two  conditions  in  Eq.  (23)  reduce  to  the  single  condition 


(25) 
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(o(t+Tj.)(t+T^)  -  a)(T)T  =  2jtn 


(26) 


f(T,Tr) 

where  existence  of  repeat  time  Tr  is  problem  dependent.  Substitution  of  Eq.  (25)  into  Eq.  (24) 
confirms  the  evaluated  frequency  function  rationality. 


a)x(x+T,.)(T+Tj.)  -  a>x(x)x  n  ^{o)(x+T^)(x+T^)  -  co(x)x}  n^ 
a)y(x+T,.)(x+Tj.)  -  a)y(x)x  ny{(o(x+T,.)(T+T,.)  -  a)(x)x}  fiy 

For  a  specific  example,  consider  m*  degree  polynomial  frequency  functions  in  Eq.  (28). 


(27) 


Wxj  =nxWi 


m  j 

cOx(x)  =  2  o)„.x 
’  i=0  *1 

m  \ 

Wv(x)  =  2  Wv  where 

y  1=0^1  Wv  =ny(0i 

m  i  yi  y  * 

0)(x)  =  2  coix 

^  i=0  ‘ 


(28) 


With  these  specific  functions,  Eq.  (26)  becomes 


Ul  I  HJ  I 

[  2  (i3j(x+Tr)  ](x+Tr)  -  [  2  WjX  ]x  =2jtn 

1=V)  l=\l 


(29) 


For  a  second  degree  polynomial  (m  =  2)  parameterized  at  zero  time  (x  =  0),  Fig.  21  shows  a  plot 
of  function  f(x,Tr)  for  a  specific  set  of  coefficient  values.  The  solution  for  the  first  (n  =  1)  and 
second  (n  =  2)  repeat  times  are  Tri  =  4.8816  and  Tri  =  7.5997.  Note  how  the  repeat  time 
increments  are  decreasing  for  each  repeat  cycle.  Corresponding  Lissajous  curve  and  individual 
x(t)  and  y(t)  responses  are  also  shown  in  Fig.  21.  Note  the  variable  frequency  Lissajous  curve  is 
closed  as  repeat  times  were  found,  but  travel  rate  at  a  given  point  along  the  curve  is  non-uniform, 
specifically  the  travel  rate  at  the  given  point  on  each  subsequent  pass  is  increasing,  even  though 
the  curve  takes  the  same  spatial  path  in  each  cycle.  Also  note  in  this  example  individual  x(t)  and 
y(t)  signals  are  not  periodic.  The  Lissajous  curve  is  spatial  periodic  but  not  temporal  periodic. 
These  observations  underscore  the  meaning  that  Eq.  (26)  is  a  repeat  condition;  not  a  periodic 
condition.  As  a  final  point  for  this  example,  Fig.  22  shows  results  for  the  case  when  cox(t)  does 
not  satisfy  Eq.  (25),  specifically,  coxi  =  0.025  0.02  =  nxcoi.  Fig.  22  shows  how  this  small 

parameter  perturbation  leads  to  distinct  functions  fx(x,Tr)  and  fy(x,Tr)  which  cross  the  27i  and  47i 
lines  at  different  points  causing  repeat  time  non-existence.  Further,  the  small  perturbation 
produces  an  open  Lissajous  curve  which  never  returns  to  pass  through  the  start  point  (x(0),y(0)) 

=  (1,0.8660)  =  (1,3  *^2/2) 


Approved  for  public  release;  distribution  is  unlimited. 

37 


Figure  22.  Variable  Polynomial  Frequency  Non-Repeating  Lissajous  Curve 
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As  another  example,  consider  sinusoid  frequency  functions  in  Eq.  (30). 


m 

ci)x(t)  +  .2  o)x,sin(iA.T) 

m  ^Xi 

WvIt)  =Wv«  +  2  cOv.sin(iXT)  where  (30) 

>0  i=i  yt  Wv  =ny(Oi  ^  ^ 

m  yi  y  ‘ 

tt>(x)  =«o  +  .2jWjSin(iX.x) 

Under  these  specific  functions,  Eq.  (26)  simplifies  to 

LcOq  +  .2  (OiSin{i?^.(x+T(.)}](T+Tj.)  -  [wq  +  .2  o)iSin{iX,T}]T  =2jtn  (31) 


For  a  single  harmonic  (m  =  1),  Fig.  23  shows  function  f(T,Tr)  for  a  specific  set  of  coefficient 
values  and  the  zero  time  parameterization  (t  =  0),  resulting  in  the  first  and  second  repeat  times 
Tri  =  6.2832  =  2k  and  Tr2  =  12.5664  =  4k.  Note  how  the  repeat  time  increments  are  constant  for 
each  repeat  cycle.  Corresponding  Lissajous  curve  and  individual  x(t)  and  y(t)  responses  are  also 
shown  in  Fig.  23.  Individual  x(t)  and  y(t)  signals  are  still  non-periodic  but  they  do 
simultaneously  repeat  their  values  at  times  x  =  0,  Tri,  Tr2,  ....  Also  note  the  variable  frequency 
Lissajous  curve  is  again  closed  and  the  instantaneous  travel  rate  at  a  given  point  along  the  curve 
is  still  non-uniform,  but  the  average  travel  rate  for  one  cycle  is  constant.  The  Lissajous  curve  is 
both  spatial  and  temporal  periodic.  To  demonstrate  that  rationally  proportional  frequency 
functions  (see  Eq.  (25))  are  not  necessary  for  a  closed  Lissajous  curve,  consider  the  results 
presented  in  Fig.  24  occurring  from  a  small  perturbation  where  cox(x)  does  not  satisfy  Eq.  (25), 
specifically,  coxi  =  0.25  4^  0.2  =  nxcoi.  Fig.  24  shows  how  this  parametric  change  causes  functions 
fx(T,Tr)  and  fy(T,Tr)  to  become  distinct,  yet  still  cross  the  27i  and  4k  lines  at  uniformly  spaced 
identical  points  facilitating  repeat  time  existence.  This  feature  gives  rise  to  interesting  behavior 
exhibited  by  the  Lissajous  curve  shown  in  Fig.  24.  The  Lissajous  curve  passes  through  the  repeat 
point  (x(0),y(0))  =  (1,0.8660)  =  (1,3 '^^/2)  every  Ax  =  2k  units  of  time  but  never  takes  the  same 
spatial  path  in  each  cycle.  Although  the  Lissajous  curve  has  remained  temporal  periodic  and 
closed  under  the  parameter  change,  the  curve  is  not  spatial  periodic.  In  other  words,  the 
fundamental  requirement  for  a  repeat  point  is  rationality  of  frequency  function  differences 
evaluated  at  two  specific  times',  not  all  times. 

Second,  consider  the  case  of  variable  phase  but  constant  amplitude-frequency,  or 

X(t)  =axCOS{C0^T-(l)x(T)} 

y(t)  =aySin{a>yT-Kt)y(T)}  ^  ^ 

Parametric  conditions  are  again  sought  where  Eq.  (22)  is  satisfied:  x,y  coordinate  values  return  to 
previous  values  after  repeat  time  Tr.  Angular  conditions  from  Eq.  (22)  after  canceling  out  equal 
frequency  terms  are 


WxT|.  -  (1)^(t-i-T|.)  h-  (Ji^fx)  =  (2jm)n^ 
WyTr  -I-  (j)y(x-l-Tj.)  -  ct)y(x)  =  (2jtn)ny 


(33) 
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Figure  23.  Variable  Sinusoid  Frequency  Repeating  Lissajous  Curve 


Figure  24.  Variable  Sinusoid  Frequency  Spatial  Non-Repeating  Lissajous  Curve 
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Repeat  time  Tr  must  again  satisfy  both  conditions  for  repeating  x,y  values.  Also,  the  closed  curve 
requirement  again  means  the  phase  function  difference  ratio  evaluated  across  the  repeat  time 
must  be  rational, 


Ilx  _  -  <l>x('^+Tr)  +  .34^ 

n  y  WyTr  +  (|)y(T+Tj.)  -  C|)y(T) 

otherwise  the  curve  is  open.  The  variable  phase  repeat  conditions  here  in  Eqs.  (33)-(34)  have 
similarities  to  the  variable  frequency  repeat  conditions  in  Eqs.  (23)-(24),  but  subtle  differences 
are  present.  Eqs.  (33)-(34)  strongly  restrict  the  class  of  time  dependent  phase  functions  that  are 
capable  of  producing  a  closed  curve,  specifically,  (|)x(t)  and  must  be  rationally  proportional 
and  out-of-phase  to  one  another.  Also  note  the  constant  frequencies  cox,  coy  are  still  present  in 
Eqs.  (33)-(34)  and  must  also  share  the  same  integer  proportionality  for  satisfying  the  closed 
curve  requirement.  This  frequency-phase  coupling  was  not  present  in  Eqs.  (23)-(24)  where 
constant  phase  angles  (t)x,  (|)y  canceled  out  from  the  expressions.  A  class  of  sufficient  time 
dependent  phase  functions  and  constant  frequencies  that  satisfy  rationality  using  base  terms  (|)(t) 
and  CO  are 


(|)^(t)  =-n^(l)(T)  =n^co 

(t)y(T)  =+  nyCt>(t)  COy  =nyCO 

In  this  case,  the  two  conditions  in  Eq.  (33)  reduce  to  the  single  condition 

(oTj.  +  <|)(t+Tj.)  -  (|)(t)  =  2jTn 


(35) 


(36) 


where  existence  of  repeat  time  Tr  is  again  problem  dependent.  Substitution  of  Eq.  (35)  into  Eq. 
(34)  confirms  the  evaluated  phase  function  rationality. 

~  ^  nx{coTr  +  (|)(T+Tr)  -  (|)(t)}  ^  n^ 

COyTr  +  (l)y(T+Tr)  -  (t)y(T)  ny{coTr  +  (l)(T+Tr)  -  (|)(t)}  Oy 


Examples  of  specific  phase  functions  are  not  considered  here  due  to  the  similarities  between  the 
variable  phase  and  variable  frequency  conditions  for  closed  Lissajous  curves. 

Third,  consider  the  case  of  variable  amplitude  but  constant  frequency-phase,  or 


x(t)  =ax(x)cos{WxT^-<})x} 
y(t)  =ay(x)sin{WyT+(l)y} 


(38) 


This  third  case  is  fundamentally  different  from  the  first  two  cases  in  that  the  time  dependent 
parameter  lies  outside  the  Lissajous  cosine  and  sine  function  arguments.  Applying  the  Eq.  (22) 
repeat  condition  for  this  varying  amplitude  case  yields 


aj^(x-HTj)cos{cOx('c+Tr)-<l)x}  =  ax(x)cos{a)xX-(j)x} 

ay(x-HTj)sin{a)y(x-HTj.)+(i)y}  =  ay(x)sin{(i)yX+<|)y} 


(39) 


For  a  specific  x,  Eq.  (39)  again  represents  two  necessary  conditions  that  a  single  value  of  Tr  must 
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simultaneously  satisfy.  Requiring  the  amplitude  and  harmonic  factors  to  repeat  independently, 
and  after  canceling  out  equal  frequency  and  phase  terms,  Eq.  (39)  leads  to  the  sufficient 
conditions 


aj^Cx+T^)  -  a^Cx)  =0 


fxCt.Tr) 

ay(x+Tr)  -  ay(T)  =0 


and 


to^Tr  =(2jrn)nx 

(OyTj.  =(2jtn)ny 


(40) 


Although  Eq.  (40)  represents /our  separate  conditions,  this  expansion  helps  to  further  develop  the 
framework.  Due  to  the  simple  nature  of  the  angular  conditions,  frequencies  cox,  Wy  must  be 
rationally  proportional,  or  integer  scalings  of  a  base  frequency  co. 


X 

3 

X 

C 

1 

X 

3 

3 

X 

c 

II 

COy 

ny 

’  COy 

3 

c 

II 

(41) 


Consequently,  repeat  time  Tr  must  be  integer  multiples  of  27i/co,  or 


The  remaining  amplitude  conditions,  with  the  known  repeat  time  from  Eq.  (42),  become 


ax(i:+Tr)  -ax(T)  =0 
ay(x+Tr)  -  ay(x)  =0 


where  Tj.  =  ^  n 


(42) 


(43) 


Note  the  harmonic  factors  have  constrained  the  amplitude  repeat  behavior  to  be  synchronized  to 
the  harmonic  frequencies  with  constant  repeat  time  increments.  Unlike  the  harmonic  factors  in 
Eq.  (39)  which  are  fundamentally  driven  by  the  same  periodic  behavior,  the  amplitude  factors  in 
general  do  not  have  any  underlying  basis,  and  their  conditions  in  Eq.  (43)  cannot  be  further 
manipulated  without  considering  specific  amplitude  functions. 

Consider  m*  degree  polynomial  amplitude  functions  in  Eq.  (44)  as  an  example. 

m 

'm  (44) 

ay(T)  +  n(T-ry.) 

In  most  situations,  these  amplitude  functions  would  never  repeat,  leading  to  an  open  Lissajous 
curve.  Exceptions  to  this  behavior  do  exist  over  a  finite  time  span,  or  indefinitely  if  piecewise 
polynomials  are  employed.  Consider  third  degree  piecewise  polynomial  amplitude  functions  that 
have  roots  at  three  consecutive  repeat  times  indicated  by  n,  n+1,  n+2  where  the  applicable  time 
span  is  Trn  <  T  <  Trn+2. 
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(45) 


ax('f)  +  ax^{T-Trn)(T-Trn+,)(T-Trn+2) 

"  ^yo  ®y  n+lX'^-^r  n+2) 

Periodic  amplitude  functions  are  easily  constructed  from  a  sequence  of  theses  polynomials.  Two- 
piece  amplitude  functions  for  ax(T)  and  ay(T)  starting  at  the  reference  repeat  time  (Tro,Tri,Tr2  and 
Tr2,Tr3,Tr4)  arc  Utilized  in  this  example.  Fig.  25  shows  a  plot  of  amplitude  difference  functions 
fx(T,Tr)  and  fy(T,Tr)  (which  are  equal  here)  parameterized  at  zero  time  (x  =  0)  for  a  specific  set  of 
coefficient  values.  Fig.  25  confirms  the  a  priori  repeat  times  (see  Eq.  (43))  whereby  amplitude 
functions  cross  the  zero  axis  at  integer  multiples  of  In/m  (Tri  =  271,  Tr2  =  471,  Tr3  =  671  where  co  = 
1).  Corresponding  Lissajous  curve  and  individual  x(t)  and  y(t)  responses  are  also  shown  in  Fig. 
25  for  approximately  3.2  revolutions.  Note  the  variable  amplitude  Lissajous  curve  is  closed  as 
repeat  times  exist.  The  curve  passes  through  the  repeat  point  (x(0),y(0))  =  (1,0.8660)  =  (l,3'^^/2) 
every  Ax  =  27i  units  of  time  but  alternates  between  two  paths  on  every  cycle  due  to  the 
polynomial  amplitude  frequency  (27r/(Tr2-Tro)  =  27i/47i  =  1/2)  being  one-half  of  the  base 
harmonic  frequency  (co  =  1).  Travel  rate  at  a  given  point  is  uniform  (fixed  frequencies),  as  can  be 
seen  from  the  periodic  x(t)  and  y(t)  signals  in  Fig.  25.  The  Lissajous  curve  is  both  spatial  and 
temporal  periodic.  As  a  final  point  for  this  example.  Fig.  26  shows  results  for  the  case  when  ax(x) 
and  ay(x)  are  perturbed  from  the  repeating  piecewise  polynomial  structure  in  Eq.  (45). 
Specifically,  when  x  >  Tr2,  additive  exponential  amplitude  perturbations  are  activated,  as 

described  below  where  ax6  =  0.005,  ayS  =  0.004,  X  =  0.5. 


=axn  +  ^x,,('f-Tr2)(''-Tr3X'C-Tr4)  +  axTe^^'^  ^r2)_l) 

T  .1  (46) 

ay(T)  =ayp -i-ay|^(T-Tj.2)(x-Tr3)(x-Tr4)  +  ay^(e^^  r2l_i) 

Amplitudes  ax(x)  and  ay(x)  do  not  satisfy  Eq.  (43)  for  x  >  Tr2.  Fig.  26  shows  how  this  exponential 
perturbation  leads  to  distinct  functions  fx(x,Tr)  and  fy(x,Tr)  after  x  exceeds  Tr2  eliminating 
existence  of  Trs  (zero  crossing  removed).  After  passing  through  the  repeat  point  at  the  end  of  the 
second  cycle,  the  perturbation  produces  an  open  Lissajous  curve  which  never  returns  through  this 
point.  The  Lissajous  curve  has  become  non-periodic  in  space  and  time. 

As  a  second  example  of  variable  amplitude,  consider  sinusoid  amplitude  functions  in  Eq. 
(47). 


m 

a„(T)  =a„  -I-  2  a,  .sin(iXT) 

0  i=l  .  (47) 

ay(x)  =ayQ  +  ,2  ay  .sin(i?^x) 

A  single  harmonic  (m  =  1)  is  utilized  in  this  example.  Fig.  27  shows  a  plot  of  amplitude 
difference  functions  fx(x,Tr)  and  fy(x,Tr)  (which  are  equal  here)  parameterized  at  zero  time  (x  =  0) 
for  a  specific  set  of  coefficient  values.  Fig.  27  confirms  the  a  priori  repeat  times  (see  Eq.  (43)) 
whereby  amplitude  functions  cross  the  zero  axis  at  integer  multiples  of  27i/co  (Tri  =  2k,  Tri  =  4n, 
Tr3  =  6k  where  co  =  1).  Corresponding  Lissajous  curve  and  individual  x(t)  and  y(t)  responses  are 
also  shown  in  Fig.  27  for  approximately  3.2  cycles.  Note  the  variable  amplitude  Lissajous  curve 
is  closed  as  repeat  times  exist.  The  curve  passes  through  repeat  point  (x(0),y(0))  =  (1,0.8660)  = 
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Figure  25.  Variable  Polynomial  Amplitude  Repeating  Lissajous  Curve 


Figure  26.  Variable  Polynomial  Amplitude  Non-Repeating  Lissajous  Curve 
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(l,3'^^/2)  every  Ax  =  2k  units  of  time  but  again  alternates  between  two  paths  on  every  cycle  due 
to  the  sinusoid  amplitude  frequency  (X  =  1/2)  being  one-half  of  the  base  harmonic  frequency  (co 
=  1).  Travel  rate  at  a  given  point  is  uniform  (fixed  frequencies),  as  can  be  seen  from  the  periodic 
x(t)  and  y(t)  signals  in  Fig.  27.  The  Lissajous  curve  is  both  spatial  and  temporal  periodic.  As  a 
final  point  for  this  example,  Fig.  28  shows  results  for  the  case  when  ax(x)  and  ay(T)  are  perturbed 
from  the  sinusoidal  structure  in  Eq.  (47).  Specifically,  when  x  >  Tr2,  multiplicative  quadratic 

secular  amplitude  perturbations  are  activated,  as  described  below  where  ax6  =  0.05,  ay6  =  0.04. 

2 

ax(T)  =axn +  ax,sin(>^T)  {1 -I- axAT-Tr2)  } 

°  *  ^2  (48) 

ay(x)  =ay^ -nay  jSin(>LT)  {1  h- ay^(T-Tj.2)  } 

Even  with  perturbation,  amplitudes  ax(x)  and  ay(x)  continue  to  satisfy  Eq.  (43)  for  x  >  Tr2.  Fig.  28 
shows  how  this  secular  perturbation  leads  to  distinct  functions  fx(x,Tr)  and  fy(x,Tr)  after  x  exceeds 
Tr2  but  existence  of  Trs  remains  (zero  crossing  retained).  After  passing  through  the  repeat  point  at 
the  end  of  the  second  revolution,  the  perturbation  produces  a  spatially  non-repeating  but  still 
closed  Lissajous  curve.  On  every  revolution,  the  Lissajous  curve  follows  a  different  path,  but  the 
curves  all  return  to  pass  through  the  repeat  point.  The  Lissajous  curve  has  become  spatial  non¬ 
periodic  but  still  temporal  periodic. 

This  section  has  explored  necessary  and/or  sufficient  conditions  for  closed  Lissajous  curves 
with  either  varying  frequency,  varying  phase,  or  varying  amplitude,  separately,  for  two 
dimensions.  Extensions  for  simultaneous  parametric  variations,  extensions  to  three  dimensions, 
and  allowance  for  multi-harmonic  terms  is  straightforward  but  requires  more  effort.  In  the 
simpler  case  of  constant  parameters,  conditions  reduce  to  the  fixed  Lissajous  curve  conditions 
noted  in  Section  4.  Fundamentally,  the  closed  curve  conditions  describe  repeating  Lissajous 
curve  points.  Analysis  has  shown  that  a  closed  curve  is  not  necessarily  a  periodic  curve  in  the 
traditional  sense  of  both  time  and  space.  By  example,  the  closed  condition  was  shown  to  expose 
temporal  periodic  curves,  spatial  periodic  curves,  and  temporal-spatial  periodic  curves. 
Additional  framework  would  be  required  to  assess  the  hierarchy  of  periodic  curve  type,  such  as 
constant  repeat  time  increments,  applicability  of  a  specific  repeat  time  across  spatially  varying 
curve  points,  or  curve  continuity  through  derivative  considerations. 
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Figure  27.  Variable  Sinusoid  Amplitude  Repeating  Lissajous  Curve 


Figure  28.  Variable  Sinusoid  Amplitude  Spatial  Non-Repeating  Lissajous  Curve 
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4.4  SECOND  ORDER  PERIODICITY 


Return  to  the  topic  of  periodicity  based  on  second  order  ROE.  For  this  development,  the 
expanded  formulation  of  second  order  ROE  employing  variable  amplitude  but  fixed  frequency 
and  phase  in  Section  5  is  chosen.  Extra  complexity  that  was  not  considered  in  Section  7,  in  the 
form  of  polynomial  center  point  terms,  multiple  harmonics,  and  three-dimensional  coordinates,  is 
present  in  the  second  order  ROE  motion  expressions.  These  features  will  be  addressed  in  a 
systematic  process.  Parametric  conditions  for  relative  motion  evolution  where  coordinate  values 
x,y,z  simultaneously  return  to  pervious  values  after  repeat  time  Tr  are  developed  from  Eq.  (49), 
which  is  an  extension  of  Eq.  (22)  from  Section  7.  Note  the  meaning  of  parametric  conditions 
here  refers  to  deputy  relative  initial  states  which  reside  in  a  layer  of  the  motion  model  below  the 
amplitude-ffequency-phase  parameters. 


x(t+Tr)  =x(t) 

y(t+T,)  =y(t)  (49) 

z(t+Tr)  =z(t) 

Using  the  expanded  structure  of  Eq.  (10),  these  necessary  conditions  become 

Xc(T:+Tr)  +  Xhi(T+Tr)  +  XshiCx+Tr)  +  Xh2(i;+Tr)  =  Xc(t)  +  Xhi(T)  +  x^hiCt)  +Xh2(T:) 

y c(t-i-Tr)  -I-  y h  1  (t+Tr)  +  y  i (r+T^)  +  y h2('c+Tr)  =  y  ^(t)  +  y  h  i  ("c)  +  y  sh  ]  ("f)  +  V  (50) 

Zc(T:+Tr)  +  Zhi(T+Tr)  +  Zshi(T:+Tr)  +  Zh2(T:+Tr)  =  Z(.(t)  +  Zhi(T)  +  z^^jiT)  +Zh2(T) 

A  set  of  sufficient  conditions,  that  require  secular  polynomial  center  point  coordinates  (xc,yc,Zc), 

first  harmonic  coordinates  (xhi,yhi,zhi),  secular  first  harmonic  coordinates  (xshi,yshi,Zshi),  and 
second  harmonic  coordinates  (xh2,yh2,zh2),  to  repeat  separately,  established  from  Eq.  (50),  is 

X  (.(XH-T,)  =  X  ^.(t)  ,  X  h  ]  (XH-Tr)  =  X  h  ]  (x)  ,  X  ^1, 1  (x+T,)  =  X  I  (x)  ,  X  h2('t+Tr)  =  x  h2(x) 

yc(x-i-Tr)  =yc(T)  ,  yhiCt+Tj)  =  yhi(x)  ,  yshiCx+T^)  =  y^hilx)  ,  yh2('f+Tr)  =yh2(t)  (51) 

Zc(x+Tr)  =  Zc(x)  ,  Zhi(x+Tr)  =  Zhi(x)  ,  Zshl('t+Tr)  =  Zshl('t)  >  Zh2('C+Tr)  =Zh2('C) 

First,  consider  the  center  point  conditions.  Progressing  to  the  ROE  secular  polynomial 

coefficient  level,  Eq.  (51)  requires 

^  x,s  1  ^  ~  ^  x,sO(  ^  ^  x,sl  ^  x,s2{(^o^)  ) 

l^y,sl{llo(''^+Ti.)}  =  b  )  +  b  y^sl  (52) 

^  z,s0(  ^ ~  ^  z,s0(  ^ ) 

By  taking  each  polynomial  term  separately,  repeat  conditions  become 
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b  x.sO'f-  ^ 

b  x,s0‘f  ^ 

=  0 

bx,sl{"o('^+Tr» 

-  bx,si{"o'f} 

=  0 

bx,s2{["o('^+Tr)]^} 

-  bx,s2{("o'^)^} 

=  0 

by,s0{U 

0 

1 

=  0 

by,sl{'io(T:+Tr)} 

-  by,s|{n„T} 

=  0 

bz,.o{l} 

0 

N 

1 

=  0 

Moving  on  to  first  and  second  harmonic  terms, 
coefficients. 


^  0=0 

-»■  No  Condition  on  b 
^  bx,si  =0  or  Tr  =  0 

^  bx,s2{"o(2'C+Tr)Tj,}  =  0 

b  j,  52  =  0  or  2t+Tj.  =  0  or  Tj.  =  0  (53) 

^0=0 

^  No  Condition  on  by  jQ 

^  by,5l{n^,T,>  =  0 

-»■  by 5|  =  0  or  Tj.  =  0 
^  0=0 

-»■  No  Condition  on  b  ^ 

.  (51)  imposes  several  conditions  on  the  ROE 


a^hlCOs{n^,(T+Tr)  -c|)^hl}  =ax,hlcos{n„x 

ay, ,|Sin{n„(T+Tr)  +(t)y,hi}  =  ay,hiSin{n,,T  +<|>y,h|} 

and  (54) 

axli2COs{2n^,(T+T,,)  -<t)x,h2}  =a^h2COs{2n„x -t|)x^2> 
ayj,2sin{2no(T+Tr)  +4>y2i2}  =3y,h2®i'^{2nQX  +<l)y^h2} 

az, h2sin{2no(x+T,)  +<l’z,h2> 


Note  these  conditions  are  the  extension  of  Eq.  (39)  in 
harmonic  factors  separately  leads  to  repeat  conditions 

Section  7.  Treating  amplitude  and 

^i,hl 

0 

II 

1 

0  =  0 

where  i  =  x,y,z 

No  Condition  on 

^i,hl 

No  Condition  on 

bi,hl 

no(x+Tr) 

and 

-  n^x  =  (2nn)l 

n^Tj  =  2jm 

(55) 

^ih2 

0 

II 

1 

0  =  0 

where  i  =  x,y,z 

No  Condition  on 

^i,h2 

No  Condition  on 

b  i,h2  ’  ^ i,h2 

2no(x+Tr) 

-  2n(jX  =  (2jin)2 

HqTj  =  2jrn 

In  Eq.  (55),  n  again  represents  a  common  integer  number  of  whole  revolutions  over  time  Tr  as  in 
Eq.  (40),  and  integers  like  nx,  ny,  nz  are  not  required  here  since  deputy  motion  along  the  LVLH 
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axes  oscillates  at  the  common  chief  mean  motion  frequency  no  or  double  frequency  2no.  Finally, 
consider  the  secular  first  harmonic  terms.  Repeat  conditions  on  the  ROE  coefficients  from  Eq. 
(51)  require 


ax,shl(T^+T,)cos{n^,(T4-Tr)  =ax,.shlWcos{n^,T 

az,shl('f+Tr)sin<n^(T+Tr)  +<t)z,,hl}  =az,shl('^)sin{noX +(1)^  shl  > 

Again  taking  the  amplitude  and  harmonic  factors  separately  gives 

('^+Tr)  -  a i,sh  1  ('^)  =  0  ^  {b I  +c 1  >  ' o^r)  =  0  where  i  =  x ,y ,z 

^  bi^,hl=0 =0  or  T,  =  0  (57) 

n^lx+Tj)  -  n^T  =  (23Tn)l  ^  n(jTj.=  23i:n 

Three  possible  cases  for  repeating  x,y,z  coordinates  exist  and  are  summarized  below. 


Trivial  Case: 

(T„n) 

=  (0,0) 

Specific  Case: 

Tr 

=  ^n 

,  T 

=  - 

-  T 

2  ■■ 

bx,sl 

=  0 

.  by, si 

=  0 

i  ,sh  1 

=  0 

>  Ci,shl 

=  0 

where  i  =  x,y,z 

General  Case: 

Tr 

=  ^n 

bx,s2 

=  0 

bx,sl 

=  0 

’  by, si 

=  0 

^i,shl 

=  0 

=  0 

where  i  =  x,y,z 

The  trivial  case  is  not  useful.  Although  the  specific  case  describes  a  correct  physical  motion 
where  the  combined  parabolic  center  point  travel  and  oscillatory  harmonic  travel  generate  a 
single  repeat  of  x,y,z  coordinates  (non-periodic)  associated  with  a  specific  time  r,  this  case  has 
limited  utility  and  is  not  of  interest  here.  The  general  case,  which  allows  for  an  infinite  number 
of  repeat  cycles  and  no  restriction  on  time,  is  the  case  of  interest  for  constructing  initial  condition 
relations  that  result  in  deputy  periodic  motion,  both  spatial  and  temporal,  under  second  order 
accuracy.  Nine  separate  conditions  are  part  of  the  general  case  listed  in  Eq.  (58)  in  terms  of  ROE 
coefficients,  which  could  have  been  alternatively  derived  from  the  second  order  compacted 
formulation  (Section  6),  or  the  original  QV  motion  expressions  (Eq.  (9)).  Only  seven  of  these 
conditions  are  potentially  independent  (see  Eq.  (11)),  which  after  simplification,  are  identified 
below  in  terms  of  deputy  initial  conditions. 
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bx^2-  (2no*0  +  yo)^  =0 

(2no*0  +  yoK^oyo  -  2xo)  =0 

by^i:  2noRo(2noXo  +  yo)  +“o(11^0  +  2yo  +  zo) 

+  (xo  +  4yo  +  zo)  +  2no(7xoyo  -  yo^o)  =0  (59) 

c x^hl  or  b  y^l :  (2n qXo  +  yo)(3n qXq  +  2y q)  =0 

Cy^hl  or  b^^i:  (2noXo  +  yo)xo  =0 

^zm’’  (2noXo  +  yo)2o  =0 

Oz^hl^  (200^0 +  yo)zo  =0 


Note  six  of  the  seven  equations  have  been  faetored,  with  one  faetor  being  the  linear  CW  drift 
condition  term  which  appears  in  all  six  factored  equations,  as  well  as  in  the  seventh  unfactored 
equation.  The  bx,s2  =  0  condition  in  Eq.  (59)  requires 


^0  =  -2iio*0 


(60) 


in  all  second  order  periodic  solutions.  Thus,  from  the  seven  conditions  in  Eq.  (59),  only  two  of 
these  conditions  are  ultimately  independent,  which  after  simplification,  are  identified  below. 


t>x,s2  »•••  : 


2noXo  +  yo  =0 

2yo  +  Zq)  (^0  '^yo  2o)  +  2nQ(7xoyo  -  yo^o) 


(61) 


In  the  CW  first  order  solution,  only  one  condition  exists  between  xo  ,70  for  a  periodic  orbit, 
which  is  the  first  condition  in  Eq.  (61).  Note  how  the  QV  second  order  solution  requires  an  extra 
condition  for  periodicity  involving  all  six  deputy  initial  conditions  in  Eq.  (61).  One  way  that  the 
by, si  =  0  condition  constrains  the  possible  periodic  motion  is  that  deputy  cross-track  or  z  axis 
motion  cannot  evolve  independently  as  in  a  first  order  model,  where  the  deputy  can  oscillate  in 
cross-track  with  zero  radial  (x)  and  zero  along-track  (y)  displacement  (assuming  collision  with 
the  chief  at  the  origin  is  somehow  avoided).  In  the  second  order  framework,  if  xo,  yo  ,^o,yo  are 
all  zero,  the  second  condition  from  Eq.  (61)  requires  zo  and  to  equal  zero  with  no  other 
physical  possibility.  Another  way  that  the  by, si  =  0  condition  constrains  the  possible  periodic 
motion  is  that  some  combinations  of  in-plane  initial  states  are  not  allowed.  Note  a  combination 
of  non-zero  initial  conditions  where  7xoyo  =  yo’^o  can  never  lead  to  a  second  order  accurate 
periodic  orbit  as  the  by,si  =  0  condition  would  become  inconsistent.  Even  though  an  extra 
constraint  on  the  initial  conditions  exists  in  second  order  accuracy,  an  infinite  number  of  periodic 
motions  is  still  possible  since  Eq.  (61)  represents  two  requirements  in  six  initial  conditions. 

Several  simple  analytical  solutions  for  the  initial  conditions  are  possible.  Consider 
eliminating  yo  in  the  second  condition  in  Eq.  (61)  in  terms  of  xo  using  Eq.  (60).  The  resulting 
single  periodic  condition  in  five  variables,  written  in  two  different  ways,  is  noted  below. 

aj(xj-2y§)  +  (2iioyo-xo)xo-(n5zj  +  zo)  =0 

or  (62) 

(njxj  -  x^  +  2no(xo  -  -  (njzj  +  z§)  =0 
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Treating  xo  as  a  specified  independent  variable,  if  yo  is  chosen  according  to  2yo^  =  xo^,  and  ^(i  is 
chosen  according  to  =  2noyo,  along  with  Yo  satisfying  Eq.  (60),  then  zo  and  must  be  zero 
and  an  in-plane  periodic  orbit  will  ensue.  Likewise,  if  is  chosen  according  to  =  no^xo^  and 
yo  is  chosen  according  to  noyo  =  ^o,  then  another  in-plane  periodic  orbit  is  created.  Further, 
suppose  xo,  yo,  zo  are  specified  and  is  chosen  according  to  =  2noyo,  along  with  yo  satisfying 
Eq.  (60),  then  Eq.  (62)  provides  a  condition  for  computing  a  corresponding  that  produces  a 
three-dimensional  out-of-plane  periodic  motion  to  second  order  accuracy.  Ultimately,  any 
specification  of  four  initial  conditions  from  the  five  member  set  xo,yo,zo,’^o,^o  with  the  fifth 
member  computed  from  Eq.  (62),  that  results  in  a  consistent  and  physically  realizable  initial 
relative  ephemeris,  leads  to  periodic  motion.  Conditions  in  Eq.  (61)  could  be  useful  for 
initializing  numerical  shooting  methods  or  differential  correction  methods  to  establish  exact 
three-dimensional  periodic  orbits  under  two-body  circular  reference  assumptions. 
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5.0  CONCLUSIONS 


An  analytical  framework  has  been  developed  for  the  eoncept  of  second  order  relative  orbital 
elements,  extending  ideas  that  first  originated  by  reformulating  the  first  order  Clohessy-Wiltshire 
Cartesian  relative  motion  solution  in  a  more  geometrieally  meaningful  characterization.  The 
framework  is  built  around  the  seeond  order  Quadratic  Volterra  relative  motion  solution,  and  is 
thus  limited  by  two-body  eircular  reference  eonditions  and  a  bounded  separation  between  deputy 
and  chief  satellites.  The  new  relative  orbital  elements  are  defined  in  such  a  way  that  geometric 
insight  to  the  relative  orbit  size,  shape,  orientation,  and  satellite  loeation  along  the  orbit  in  three- 
dimensional  settings  is  preserved.  The  new  element  set  is  based  on  amplitude  and  angular 
quantities  related  to  the  eoncept  of  an  instantaneous  ellipse.  Through  the  use  of  auxiliary  cireles, 
a  given  element  set  can  be  used  to  geometrieally  eharacterize  the  satellite  orbit.  Seeond  order 
motion  produees  varying  amplitude  and  phase  behavior,  and  the  new  relative  orbital  element 
variables  are  well-suited  for  deseribing  this  complex  behavior.  The  new  framework  is  highly 
related  to  the  coneept  of  variable  Lissajous  space  curves,  which  was  utilized  to  develop  the  new 
element  variables.  New  relations  were  developed  to  deseribe  and  quantify  how  first  order  relative 
orbital  elements  are  distorted  by  seeond  order  gravitational  effects,  that  may  be  useful  in 
dynamieal  guidance  and  estimation  error  analysis.  Finally,  new  eonditions  for  periodic  orbits  to 
second  order  aeeuraey  were  developed  from  the  new  element  set  and  repeating  Lissajous  eurve 
conditions.  These  seeond  order  eonditions  were  shown  to  be  more  restrietive  on  initial  condition 
sets  that  produee  periodic  motion,  when  compared  to  eorresponding  first  order  eonditions. 
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6.0  RECOMMENDATIONS 


Several  follow-on  researeh  and  development  efforts  to  the  work  presented  in  this  report 
should  be  pursued.  First,  a  software  package  that  utilizes  computer  generated  three-dimensional 
color  visualizations  to  facilitate  the  practical  use  of  the  second  order  relative  orbital  elements  and 
the  geometrical  information  they  depict  should  be  developed.  A  tool  such  as  this  could  enhance 
flight  control  operations,  space  situational  awareness,  and  rapid  comprehension  and  adaptation  of 
complex  multi-satellite  formations.  Second,  the  analytical  framework  and  underlying  numerical 
calculations  of  second  order  relative  orbital  elements  should  be  applied  to  maneuver  planning 
and  execution  strategies.  Orbital  maneuver  strategy  based  on  the  new  motion  description  could 
be  exploited  to  provide  more  efficient,  more  accurate,  and  more  transparent  guidance  and  control 
functions  involving  multiple  close-proximity  cooperative  satellites.  Third,  the  new  perspective 
on  necessary  conditions  for  periodic  relative  orbits  should  be  further  examined  and  potentially 
utilized  in  orbital  flight  control  procedures  addressing  the  objectives  of  maneuvering  to  and 
from,  and  maintaining,  periodic  relative  orbits  in  an  optimal  fashion.  Increased  knowledge  from 
this  subject  area  should  also  be  considered  for  initialization  of  numerical  search  algorithms 
attempting  to  achieve  high  precision  conditions  for  achieving  periodic  relative  orbits. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


cw 

Clohessy-Wiltshire 

ECI 

Earth-Centered  Inertial 

LVLH 

Local- Vertical  Local-Horizontal 

NLS 

Nonlinear  Simulation 

OED 

Orbital  Element  Differences 

QV 

Quadratic  Volterra 

ROE 

Relative  Orbital  Elements 
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